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THE EFFECT OF STERILISED MALES ON A NATURAL 
TSETSE FLY POPULATION 


H. R. Simpson 


Rothamsted Experimental Station 
Harpenden, Herts., England 


1. Introduction 


The control of an isolated insect population by the release of arti- 
ficially sterilised males has been successfully carried out by the En- 
tomology Research Branch of the United States Department of Agri- 
culture. The apparent eradication of the screw-worm [Callitroga 
hominivorax (Cqrl.)] from the island of Curagao was described by 
Baumhover et al. [1955], and Knipling [1955], in a theoretical analysis 
of the principle, has examined the biological and economic conditions 
necessary for its success. The biological prerequisites are: 

(i) The method of sterilisation must not adversely affect the 
mating behaviour of the males. 

(ii) Adequate dispersion of the sterilised males amongst the natural 
female population must be obtainable. 

(iii) If the females mate more than once, the sperm from sterilised 
males must be capable of competing with that from normal 
males. 

Since the tsetse female normally only mates once, the question was 
raised’ of the effect of this method on a natural tsetse fly population.” 
The purpose of this paper is to present a theoretical examination of the 
method, which though based on numerical data derived from the 
observed life cycle of the tsetse fly, is sufficiently general to indicate its 
application to other insect populations satisfying the above conditions. 

In Section 2 the necessary mathematical background is sketched. 
A mathematical model of the natural tsetse population is developed 
(Section 3) and the theoretical effect of sterilised males is considered 


1By Mr. J. K. Chorley (Director of Tsetse Operations in Southern Rhodesia), through the Co- 
lonial Office Tsetse Fly and Trypanosomiasis Committee. 

2Previous experimental work on a similar idea has been reported by Vanderplank [1947]. In this 
case sterilisation of the female was effected by mating with males of a closely related species; an attempt 
to apply it in field control failed, possibly due to natural segregation of the species. In the Curacao 
experiment the insects were sterilised by gamma-irradiation and this cause of failure would not apply. 
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(Section 4). Some notes on the numerical calculations of the results 
are given in Section 5. These results are discussed fully in Section 6, 
and some general comments appear in Section 7. 


2. Theoretical Background 


The problem was treated mathematically by a technique due to 

P. H. Leslie [1945] in which the expected numbers of females in each age 

group are considered in detail. If suffices r, n refer to age and time 

respectively where the units are the same and the age origin is the 

pupation of the larva, let 

v,.,_ be the expected number of females in the age group (r, r + 1) 
at time n; 

m,., be the probability that a female in the age group (r, r + 1) alive 
at time 7 is still alive (in the age group (r + 1, r + 2)) at time 
(n + 1); 

¢,,, be the expected number of female pupae deposited in the time 
interval (n, n + 1) by a female in the age group (r, r + 1) which 
are alive at time (n + 1) [this definition of ¢,,, enables the first 
equation of formulae (1) to be simply stated, but it must be taken 
into account when determining the survival factor p, defined 


below]: 
then 
Vo,n+1 
(1) 
= Tr in (r => 


In matrix form these equations may be written 
= 
where v, is the column vector and 


If P, = P is independent of n then v,..,—pv, where p is the latent root 
of greatest modulus of P, and it can be shown that for matrices of this 


‘form (if the elements are such as to be biologically meaningful) p is real 


and positive. .Thus the population level tends to a constant value only 
ifp = 1. 
In order to make this so for the model here considered, the elements 
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of P need to be adjusted accordingly (see Sec. 3). Strictly speaking, 
the population equilibrium is theoretically neutral for p = 1 but unstable 
for p differing from unity by an arbitrary small quantity. Genuine stable 
equilibria in mathematical models ef population systems only appear 
possible where density dependent assumptions are incorporated into the 
system, but this would be a refinement beyond the intention of the 
present study. 

Biologically the assumption p = 1 may be justified if it is supposed 
that under the conditions prevailing the effects of birth and death are 
so closely equated that for a reasonable time interval no marked popula- 
tion change need be expected. If we are prepared to limit our con- 
siderations to a reasonably short finite period of time, the population 
may be regarded as virtually constant. 

The introduction for a limited period of disadvantageous conditions 
(reflected in new values of ¢ or x) would depress p below unity and result 
in a population decline, so that on the cessation of those disadvantageous 
conditions and the return of p to unity, the population would assume a 
new and lower ‘stationary’ level. ‘The particular disadvantageous 
condition considered in this paper is the temporary reduction in birth 
rate brought about by the introduction for a limited period of time of 
sterilised males. 


3. The Natural Population 


The development of the mathematical theory requires certain specific 
assumptions concerning the biology of the tsetse fly. Those given 
below are based on Glossina palpalis (R.D.) and G. morsitans (Westw.) 
but the variations between the different species are thought to be 
relatively unimportant. 

Assumption 1. 

The natural population is ‘stationary,’ i.e. the v, 7, @ of (1) are inde- 
pendent of n and the population numbers remain constant. In Sec. 
6.5, however, a finitely oscillating population model is considered. 
Assumption 2. 

The death process is Poisson, i.e. the probability that an individual 
dies in a small time interval (¢, t + 6¢) is proportional to the length of 
the interval, 6¢. Some such assumption is necessary in order to deter- 
mine the z,,, ; the Poisson death process agrees reasonably well with 
observation and is mathematically the simplest. 

Assumption 3. 
(i) The pupal period is 4 weeks. 

(ii) An adult male is sexually mature one week after emergence 

from the pupa and remains potent until death. 
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(iii) An adult female mates within 3} days of emergence and does 
not mate again. 

(iv) A fecund female produces one pupa 17} days after mating and 
subsequently one after every 10} days until the maximum 
reproductive age of 30 weeks (measured from emergence) is 
reached. Male and female pupae are equiprobable. 

The particular periods of time specified were chosen as convenient 
multiples of a half-week; the reproductive age limit imposed is arbitrary. 
The restriction of single mating is relaxed in Sec. 6.4. 

In what follows the time unit is one week. Let M, be the expected 
number of adult males at time n, and F,, be the expected number of adult 
females below the maximum reproductive age at time n, ie. 5°,» 45 - 
Let x, » be the male, female death rates respectively; since the death 
process is Poisson the expected durations of life T are 1/« and 1/u. 

Only effective pupae, i.e. those pupae that will eventually emerge 
as adults, are considered. A survival factor p, is therefore introduced, 
so that 


@) 
0 otherwise. (by Assumption 3). 
Moreover 
= (by Assumption 2). (3) 
oe (r > 4). 


T is then the expected duration of adult life. 

For convenience all population numbers are expressed as proportions 
of the total number of adult females in the natural population; thus, 
omitting the suffix n, >>°., », =1. But from (1) and (3) 


{" O<r<4) (4) 
whence 
and F=1-—e™, (5) 
But from (1), (2) and (4) 
11 10 
r=0 r=2 r=2 
hence 
21 — 
p= (6) 


+e “(1 — 
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Equations (4), (5) determine the natural female age distribution; for 
the male population 


Mas eM, + V3 in (7) 
and in stability 
M =—3— = 8 
1-—e“ 1-—e"“ ®) 


4. The Effect of Sterilised Males 


It is now assumed that: 
Assumption 4. 
A constant number S of sterilised males is introduced weekly for 
a given number of weeks ¢ beginning at n = 0. 
Assumption 5. 

The mating behaviour of the male is unaffected by the sterilisation 
procedure. 
Assumption 6. 

No density dependent factors act upon the population as it declines 
(i.e. the survival factor p remains constant). 

The probability f, that a mating in (n — 3, n + }) is fertile is no 
longer unity and following Baumhover et al. [1955] we assume that: 
Assumption 7. 

f, is equal to the ratio of mature to mature plus sterilised males at 
time n. 

The original population will initially decrease; but after n = ¢ the 
number of sterilised males will decline exponentially until their effect is 
negligible and the original age distribution will eventually reappear at a 
lower level. A measure of the effect of the method is therefore the 
proportion of the original population to which the disturbed population 
falls. 

Let S, be the expected number of sterilised males at time n, and 
«’ be the sterilised male death rate. 


Then 
(M,, ¥s.0-1) 
where 
So = 0 
and 


(10) 


Mati ™ 
(n > 0) 
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(These equations may be solved explicitly, giving 


1-e™ 
O<n<d 
(n > t) 


with a maximum at n = t of 


This may be used to determine when the sterilised male population 
becomes negligible). 
Finally 


Orin ° 


} Pod v,., iS the total number of adult females, both fertile and 
sterile, below the maximum reproductive age at time n. The number 
of fertile females at time n is given by 


Fast e"(F, = fn-29¥33,n) + (n > 0) (11) 
5. Calculation of the Results 


The calculations required to determine the reduction in the natural 
population caused by the release of sterilised males would be impossible 
without the aid of an electronic computer. The numerical results 
presented in this paper were computed on the Elliott 401 at Rothamsted, 
a description of which has been given by Lipton [1955]. 

The basic programme written for the 401 was based on the assump- 
tions set out in Sec. 3 and Sec. 4. It requires as initial data yp, x, x’, S 
(expressed as fm) and ¢t. The original population figures are calculated 
[from (4), (5) and (8)]; and from these values and the recurrence relations 
[(1), (7), (10), and (11)] the population figures for succeeding weeks are 
derived. At intervals of four weeks n, M, , von, F, and >>%%, »,,. (ie. 
time; adult male population; female pupae deposited; fecund and total 
adult female populations) are printed out, and the process is allowed to 
continue until the final level is reached (F, = )/°2, »,,.). (The com- 
puting time is approximately 12 seconds, plus a further 8 seconds for 
output, for each four weeks.) Certain modifications were later made to 
this basic programme to enable various alternative assumptions to be 
considered; these are discussed in Sec. 6. 


This machine has a word length of 32 binary digits, i.e. about 9 decimal places, and a nominal 
range of [—1, 1); the restriction eel 4¥r = 1 ensured that most of the population figures occurring 
in the computations lay within this range. The exception, S, , was expressed as 2*.Sn (0 < Sn < 1). 
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In practice the effect of round-off errors and of random machine 
errors had to be examined. It was found that round-off errors could be 
ignored, and that with regard to machine errors the process was self- 
correcting in the sense that small errors [0(10~*)] in an intermediate 
stage of the calculations were attenuated in the subsequent calculations 
and led to negligible differences in the final results. Large machine 
errors usually resulted in the appearance of a negative population value 
[due to exceeding the range (—1, 1)], when the calculations were 
restarted ab initio. 


6. Examination of the Results 


The results obtained by the computations outlined above are pre- 
sented in Tables 1-5. In Tables 1-4 the effects of the introduction of 
sterilised males are shown in terms of the expected changes in the numbers 
of females. The figures in the bodies of the tables give the numbers of 
females, at various times, expressed as percentages of the initial numbers 
of females, the time in weeks that has elapsed after the first introduction 
being represented by the symbol n. The number of sterilised males 
introduced in each of ¢ weeks is expressed as a fraction f of the initial 
number of males. 

The figures given in the tables are rounded off to two decimal places, 
except for values < 0.01 where the first significant figure is given. 
Stable levels, where established by computation, are shown in italics; 
a * indicates apparent elimination, i.e. a value of < 0.0°5. 

The results depend on the assumptions of expectation of life of the 
various types of fly. In Tables 1, 3, 4, and 5 the standard assumptions 
are made, namely that the expectation of life of females is six weeks and 
that of males, both sterilised and unsterilised, is three weeks. In Table 2 
the effects of various changes in these assumptions are shown. 


6.1 Table 1 shows the results obtained for the standard case when the 
number of sterilised males introduced (determined by f and ¢) is varied. 

The values of f that are of practical interest are limited by two 
considerations: 


(i) A substantial reduction (to, say, <1%) must be obtained within a 
reasonable time. 
For f = 0.01, 0.02, 0.03, 0.04, 0.05 and t = 96 the final levels 
are \ = 63.92, 34.41, 13.87, 3.50 and 0.46 respectively. Thus only 
f = 0.05 need be considered. 

(ii) It must be economically feasible to produce the required number of 
sterilised males weekly. The upper limit taken is f = 2. The 
results for f = 0.05, 0.1, 0.2 (¢ = 32, 64, 96) and f = 0.5, 1, 2, 
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TABLE 1 


Errect OF WEEKLY INTRODUCTION OF STERILISED MALES 
(Standard Assumptions of Expectations of Life) 


f = 06 f=. f=2 
t 32 64 96 32 64 96 32 64 96 
n 
40 51.68 42.95 42.95 26.61 18.08 18.08 8.01 4.55 4.55 
60 51.34 20.54 20.54 26.07 2.69 2.69 7.33 0.16 0.16 
80 61.85 15.11 5.95 26.08 0.58 0.10 7.384 .0°2 .0*1 
100 15.07 0.73 0.55 .0°7 07 * 
120 0.46 . 
* 
f= 6 f=2 
t 16 32 48 16 32 48 16 32 48 
n 
40 14.46 0.88 0.63 5.42 0.20 0.18 1.69 0.06 0.06 
60 0.42 .076 §.40 0.02 .0% | 1.66 .097 .0% 
80 081 0.01 .0°1 * 
100 043 
Final Levels 
t 16 32 48 64 96 
f 
.05 51.35 15.07 0.46 
af 26.08 0.55 
2 7.34 
5 14.46 0.42 .0% 
1.0 5.40 0.01 * 
2.0 1.66 


(¢ = 16, 32, 48) are shown (the values of ¢ are chosen to give a 
significant final level in the majoiity of cases); the final levels 
for these values, and others computed by the 401, are shown 
graphically in Fig. 1. 

The effect of the variations in f and ¢ is complex; the decline in the 
final levels if either is constant is approximately exponential. It is 
worth noting that if the total number of sterilised males introduced, 
ftM, is fixed, it is better to introduce a smaller number weekly for a 
longer period: the final level attained when f = 0.05, t = 96, is 0.46; 
when f = 2, ¢ = 16, the final level is 1.66, although the total number of 
sterilised males required is nearly seven times greater in the latter case. 
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figure 1 
Final A Obfa. 
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6.2 The effect of variations in the expectations of life, keeping the 
number of sterilised males introduced fixed, is shown in Table 2. The 
variations considered are: 


(i) The expectation of life of the sterilised males halved, to allow for 
a possible harmful effect of irradiation on longevity; 
(ii) The expectation of life of the female doubled, and 
(iii) The expectation of life of all flies doubled, to allow for inaccuracies 
in their estimates. 


The results for these cases, together with those for the standard case, 
for f = land t = 16, 24, 32, are given. [rom these figures the following 
conclusions may be drawn: 


(i) A relative decrease in the expectation of life of the sterilised male 
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or a relative increase in the expectation of life of the female will 
substantially reduce the efficacy of the method. 

(ii) A uniform increase in the expectation of life (preserving the 2:1:1 
ratio) leads to better results (effectively, more sterilised males are 
introduced per lifetime). Conversely, a uniform decrease will lead 
to worse results. 


TABLE 2 
EFFECT OF CHANGES IN ASSUMED EXPECTATIONS OF LIFE 


(Figures in Brackets Indicate the Expectation of Life of Females, 
Males, and Sterilised Males Respectively ) 


fel 
Standard case: (6; 3, 3) Life of sterilised males halved: 
(6; 3, 3/2) 
i 16 24 32 16 24 . 32 
n 
40 5.42 0.69 0.20 17.19 4.89 0.97 
60 5.40 0.54 0.02 17.21 4.87 0.91 
80 0.01 
Life of females doubled: Life of all flies doubled: 
(12; 3, 8) (12; 6, 6) 
t 16 24 32 16 24 32 
n 
40 13.68 "a ef 1.49 4.04 1.53 1.09 
60 13.84 3.47 0.58 2.58 0.30 0.06 
80 13.85 3.46 2.45 0.15 06 
100 2.44 13 -071 
120 .0°6 


6.3 Since the number of sterilised males introduced per week will be 
severely limited by practical considerations, some method of initially 
reducing the population numbers is desirable. We consider here a single 
application of an insecticide followed immediately by the introduction of 
sterilised males; it is assumed that the adult population is reduced to 
k% of its original size, and that S = kM/100 sterilised males are intro- 
duced per week. 

The results for k = 5, 10, 20 (under standard assumptions concerning 
the expectations of life) for t = 16, 24, 32, are shown in Table 3. For 
comparison the results for the corresponding cases in which no initial 
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reduction is attempted are given in brackets; and the first column (S = 0) 
shows the effect of the insecticide alone. 

It is clear from these figures that a single application of an insecticide 
which attacks only adult flies has comparatively little effect, owing to 
the large number of flies in the pupal state [for the case considered 
(u = $) the ratio of (effective) female pupae to adult females is roughly 
3:5]. In practice of course an insecticide is applied several times over 
a period of weeks. 


: TABLE 3 
EFFEcT OF AN INITIAL REDUCTION OF THE ADULT PoPULATION TO k% 
OF THE ORIGINAL VALUES 


(Values When There Is No Initial Reduction Are Shown in Brackets) 
The Values of f Refer to the Original Male Population 


k=6 f = 05 
n S=0 t= 16 24 32 
40 36 .06 13.49 (74.07) 7.82 (62.49) 3.93 (51.68) 
60 13.52 (74.07) 7.81 (62.44) 3.74 (51.34) 
k = 10 f= 
n S=0 t= 16 24 32 
40 89.42 8.99 (56.98) 3.56 (40.35) 1.14 (26.61) 
60 9.02 (56.99) 3.54 (40.27) 0.91 (26.07) 
3.53 (40.27) 
k = 20 f=2 
n S= t= 16 24 32 
40 46 15 6.01 (36.66) 1.538 (18.71) 0.38 ( 8.01) 
60 6.02 (36.67) 1.48 (18.60) 0.18 ( 7.33) 
80 (E34) 


The effect of the insecticide application on the results obtained from 
sterilised males varies considerably from case to case. For k = 5, 
t = 16, the ratio of the corresponding final levels is 5.48; for k = 20, 
t = 32, it is 43.2. 


6.4 The restriction [Assumption 3(iii)] that the female should mate only 
once, though in accordance with accepted ideas, has recently been 
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questioned (Nash [1955]); it is therefore desirable to examine the effect 
of relaxing this condition. 

Let (p, q) refer to the case where up to p matings are possible at 
intervals of q weeks and let F, S denote fertile, sterile matings respec- 
tively. If it is assumed that no further mating occurs after a fertile 
mating, the possible cases are, for p = 2: F; SF or SS; 

and for p = 3: F; SF; SSF or SSS, 
where, for example, SSF means two sterile matings and one fertile 
mating in that order (at intervals of g weeks). Equations (2) and (11) 
no longer hold; the necessary modifications are straightforward but will 
not be given here. 


TABLE 4 
Errect or Muttiete Marines (p Matinas at INTERVALS or g WEEKS) f = 1 
p=2,q=1 p=2,q=2 
t 16 4 32 16 4 32 
n 
40 14.64 2.46 0.41 13.33 2.07 0.34 
60 15.94 3.68 0.32 14.61 3.13 0.24 
80 15.99 3.77 0.42 14.67 3.23 0.32 
100 3.78 0.438 
p=8,q=1 p=3,q=2 
t 16 24 32 16 32 
n 
40 25.56 9.22 2.38 21.93 6.99 1.61 
60 25.1 8.59 1.80 21.50 6.52 1.14 
80 25.12 8.55 1.75 21.49 6.48 1.09 
100 1.74 
Final Levels 
t 16 24 32 
Pd 
1,1 5.40 0.54 0.01 
2,1 15.99 3.78 0.43 
2,2 14.67 3.23 0.32 
3,1 25.12 8.55 1.74 
3,2 21.49 6.48 1.09 
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Table 4 gives the results for (p, gq) = (2, 1), (2, 2), (3, 1), (3, 2), for 
the standard case. The summary table compares the final values with 
those obtained under the original assumptions [i.e. (p, g) = (1, 1)]. In 
terms of ratios, the reduction in the efficacy of the method is consider- 
ably greater for larger values of t. 


6.5 The possibility of natural (non-random) fluctuations in the normal 
tsetse population (e.g. due to seasonal climatic changes) has hitherto 
been ignored (Assumption 1). In this section a finitely oscillating 
population model, with a period of 32 weeks, will be set up. 

With » = 3, x = 3, the survival factor p is found to be 0.6453. 
If we define 


(l<n< 16) 
0.378 (17 <n < 32) 


and put p, = 1.000 and 0.378 in alternating periods of 16 weeks, then 
the population is found to oscillate; the values of F,, (the fecund female 
population) for n = 0(4)32 are given in the first part of Table 5. (These 
values of p, were found by trial and error; in practice, after setting 
up the initial stable population values, the process was allowed to run 
until n = 112 and the population values at that time were taken as 
the new initial values.) 


TABLE 5 
NATURALLY OSCILLATING POPULATION 
Final level obtained when S = M sterilised 
The oscillating population: adult males are introduced for ¢ weeks, 
female population at time n beginning at n = no 

n F, n, |t: 16 24 32 

0 1.02 0 8.97 0.94 0.05 

4 1.16 4 6.05 0.63 0.03 

8 1.29 8 4.90 0.57 0.02 

12 1.24 12 3.86 0.48 0.01 

16 1.13 16 3.49 0.38 0.007 

20 99 20 5.49 0.62 0.01 

24 90 24 6.65 0.65 0.01 
28 .94 28 8.21 0.77 0.03 

32 1.02 32 9.10 0.96 0.05 


The magnitude of the oscillation, though small, is sufficiently large 
for the present purpose which is to examine the effect of introducing 
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S = M, sterilised males weekly for ¢ = 16, 24, 32 weeks beginning at 
different points (n = mp») of the cycle. The results for no = 0(4)32 
appear in Table 5; the final levels obtained may be compared with 
those for the standard case (A = 5.40, 0.54, and 0.01 for t = 16, 24 and 
32). The optimum result is obtained when n, = 16. 


7. Discussion 


The following points must be emphasised: 


(i) The practical difficulties involved in the application of the method 
(c.f. Knipling [1955]) are outside the scope of this paper; they are 
undoubtedly very great and may well prove prohibitive. 

(ii) The assumption (6) that no density dependent factors operate is 
important; although none are known for the tsetse fly, their non- 
existence is by no means certain. 

(iii) The mathematical treatment has been entirely deterministic, i.e. 
no account has been taken of random fluctuations in the population 
numbers. 


Nevertheless, we may reasonably conclude from the results quoted 
that it is impractical to attempt the control of a high density tsetse 
population (c.1000 males per square mile) by the irradiation method, 
since a large number of sterilised males are required. For example, to 
obtain a reduction of such a population to 0.001% of its original size: 
with f = 0.1, 100 male pupae per square mile must be collected, irradi- 
ated, and released weekly for 96 weeks; with f = 2 the requirements are 
2000 males per square mile for 32 weeks. Moreover, these estimates do 
not allow for mortality amongst the collected pupae. 

Where, however, the tsetse population is of low density (c.100 males 
per square mile) the irradiation method may prove effective, especially 
when the low density is produced by an insecticidal application immedi- 
ately before the release of sterilised males. 

Experimental work on the sterilisation of tsetse flies by irradiation 
with gamma rays, now nearly completed, shows that both sexes are 
sterilised, the females from the irradiated pupae showing no ovarial 
development, and the males producing spermatazoa which fail to 
fertilise even normal females. The effect is not absolute, and the extent 
of the exceptions to it remains to be fully worked out, as does also the 
degree to which females will mate more than once. When these are 
known, it will be possible to determine the extent to which the decline 
of population resulting from the introduction of irradiated males will be 
affected by these exceptions to the assumptions on which the calculations 
have been based. 
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8. Summary 


The use of artificially sterilised males has been suggested as a means 
of controlling tsetse populations. A mathematical model of a natural 
tsetse fly population is set up, and the theoretical effect of the introduc- 
tion of sterilised males is examined. Numerical results obtained on 
an electronic computer are presented and discussed. 
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MAXIMUM LIKELIHOOD ESTIMATION FROM 
INCOMPLETE DATA 


H. O. 
Iowa State College, Ames, Iowa, U.S.A. 


I. Introduction 
If a random sample of size N has been drawn from a population 
involving unknown parameters 6, --- 6 , the latter may be estimated 


from the sample by the well-known technique of maximum likelihood 
the properties of which have been frequently examined in the past. 
The equations resulting from this technique are often simple but 
sometimes require iterative procedures which are usually supported 
by special aid-tables. 

However, situations frequently arise in which the data are ‘in- 
complete’ in the sense that the information contained in the random 
sample of size N is not completely available for estimation. The two 
best known examples of such ‘incomplete’ information have become 
known under the names of ‘truncated’ and ‘censored’ samples but there 
are numerous other situations. The literature of the past decades is 
abundant with papers dealing with special cases of such ‘incomplete 
information’ for special distributions (see e.g. Finney, D. J. [1949]). 
Many authors, considering the computational work involved in maxi- 
mum likelihood estimation prohibitive, have suggested alternative 
methods. (See e.g. Moore [1952].) 

The purpose of the present paper is to sirnplify and unify the maxi- 
mum likelihood computations of estimates from ‘incomplete data.’ 
It should be stressed that we are using here standard maximum likeli- 
hood estimation as introduced by Fisher and as used by many statis- 
ticians, and further, that we are not concerned here with discussing 
the small sample efficiency of these methods. We do, however, offer 
a new approach to the computational evaluation of these estimators 
from data in incomplete samples which represents a considerable 
improvement on the methods in current use. We shall in fact show 
how maximum likelihood estimation from ‘incomplete data’ can be 
simply reduced to that from a complete sample. The method is akin 
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to the ‘missing plot technique’ in analysis of variance, in fact the latter 
may be regarded as a special case of our simple iterative method de- 
scribed below. This procedure is applicable to any situation in which a 
maximum likelihood procedure for the complete sample is available. It 
covers a wide class of problems in which incomplete data arise incor- 
porating as special cases problems of truncation or censoring at either 
end of the distribution or in the central section or any combination 
thereof. It is not claimed that in every particular problem of truncation 
or censorship the present method is simpler than an already existing 
method specifically developed for that situation and supported by 
special aid tables, nevertheless this will be so in many cases. The 
main purpose of the paper is to provide a method of complete generality 
depending only on tables of the probabilities of the distribution and 
not on aids or methods which vary with the nature of the incompleteness 
of the data, and more particularly provide a workable method for the 
multitude of cases in which no special methods are available. We 
commence by describing the method for discrete distributions leaving 
the discussion of continuous distributions to a second paper in which 
we shall also discuss the fit of a truncated Gamma distribution to 
rainfall data which problem gave rise to the present study. 


2. Discrete distributions with missing frequencies (truncation) 


We first describe the method in terms of a simple numerical example 
which we employ to introduce the notation; we then give a general 
definition of this case and proof of the formulas used. 


Example 1. 


Snedecor [1956, p. 483] quotes data of Leggatt [1933] on the pollution . 
of seeds of Phleum pratense by the presence of a few noxious weed 
seeds. Table 1 shows the frequency with which 2, 3, 4, --- of these 
weed seeds were found in 78 quarter-ounce samples. Actually Snedecor 
also gives the frequencies of quarter ounces with 0 or 1 weed seed but 
we assume, for purposes of illustration, that these frequencies are not 
available, that we must estimate the mean number of weed seeds per 
quarter ounce from the remaining frequencies and that this count x 
follows a Poisson a distribution e~°@*/xz! We introduce the following 
notation: 


n, = set of observed frequencies >>; n; = 7, (1) 


«n; = c-th estimate of j-th missing frequency >>; .n; = .n’, (2) 
6) 


probability of z = 7 (here e~’6*/i!), (3) 


| = 
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f(j, 0) = probability of x = j (here e~’6’/j!), and (4) 
f(0) = 4). (5) 


The iterative procedure described below and illustrated numerically 
in Table 1 yields the maximum likelihood estimate of the Poisson mean @ 
(as will be shown below) and consists of the following steps: 

(ao) Estimate initial values of the missing frequencies on; (in Table 1 
oM = 4, on, = 14, so that on’ = 18). 

(bo) Using the n,; and the on; compute an initial estimate ,6 of the 
parameter from the maximum likelihood equation appropriate to 
the complete sample. In the example this equation is the mean of 
the complete distribution so that 

16 = (Doing + + on’) = (279 + 14)/96 (6) 

= 3.052. 


(a,) Using ,@ compute ‘improved’ estimates, in; , of the ‘missing 
frequencies’ by proportional allocation based on the Poisson distri- 
bution, i.e. 


nm; = nf(j, — f). (7) 
In the example (see Table 1) we have 
No = 4.564, sm, = 13.935, so that yn’ = 18.499. 


The reason for the high decimal accuracy in these preliminary 
calculations will be apparent in Section 5. 

(b,) Using the n; and the ‘improved’ ,n; , compute an ‘improved’ 
estimate .@ of @ from the maximum likelihood equation. In the 
example 


20 = (Do ins + jm)/(n + wm’) 


= (279 + 13.935)/(78 + 18.499) = 3.036. 


This value is already close to the previous one of ,@ = 3.052 but for 
higher accuracy two more cycles, shown in Table 1, yield ,@ = 3.026 
which is the maximum likelihood estimate of @ to almost 3 decimal 
accuracy. 

The convergence of the process has been found to be extremely 
rapid in some 30 examples worked, provided care is taken in the initial 
choice of the on; and provided that the truncation is not ‘too drastic,’ 
i.e. provided n’/(n + n’) is below 4. Above this value the convergence 
may still be serviceable if accelerated by such methods as described in 
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TABLE 1 
Maximum LIKELIHOOD FIT TO PoIssON DISTRIBUTION WITH 
MIssING FREQUENCIES 


Count 
t= | 
t j ni 1Nj nj | f(9, 18) 28) 3) 
v | 4 4.564 4.655 4.690 | .0473 .0481 .0484 
1 14 13.985 14.119 14.205 | .1444 .1459 .1466 
4 26 From table of Poisson 
distribution 
3 16 
4 18 
5 9 
6 3 
7 5 
8 0 
9 1 
10 | 
Totals n=78 cn’=18 18.499 18.774 18.895 |f=.1917 .1940 .1950 
ing =279| jenj=14 13.935 14.119 14.205 
e8=3.052 3.036 . 3.029 


Section 4. [In the example n’ = 19, n = 78, n’/(n + n’) = 0.2.] The 
theoretical background of the convergence will be discussed in a second 
paper and the variance computation of the resulting maximum likeli- 
hood estimates in Section 5. The effect on convergence of choosing 
poorer starting values is discussed in Section 4. 

The above is an example of what is known as a sample from a 
‘truncated distribution’ or more precisely we are speaking of a distri- 
bution truncated at the lower tail end. In such a situation we omit 
from the original distribution (here the Poisson distribution) all fre- 
quencies for which x < 7* where 7* is the point of truncation (here 
i* = 3). We then proceed to draw a random sample of fixed size n 
from the ‘truncated distribution’ consisting of the (Poisson) frequencies 
for which x > 7* divided by their total. An alternative description of 
the above process of sampling is to draw random samples of unknown 
size from the complete distribution and then to consider the subpopu- 
lation of samples for which the number of x values with x > 7* is a 
constant sample size n, whilst the remaining values with x < 7* are 
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unknown. The first description of the sampling procedure provides a 
clear theoretical background for the method of estimation of the param- 
eters of the truncated, and hence of the complete distribution (e.g. by 
maximum likelihood). The second description is supposed to explain 
how such ‘truncated samples’ may arise in practice and has recently 
been critically discussed by F. N. David and N. L. Johnson [1952], 
who raise the question whether there may be practical situations where 
this method of sampling does in fact not apply. Nevertheless it has 
been almost universally accepted in the past. 

The above case of ‘truncation’ can be described briefly as a sample 
in which all frequencies for counts x < 7* are ‘missing.’ 

In the present treatment we deal with a more general situation of 
missing frequencies: We assume that the set of all possible counts 
x = 0, 1, 2, --- is subdivided into two subsets, namely: 

(a) the counts x = 7 for which frequencies n; were observed 
(b) the counts x = j for which frequencies n; were not observed. 
(In Table 1, 7 = 0, 1 andz = 2, 3, 4, ---.) 

Again, it is assumed that >>; n; = n is fixed in repeated sampling 
whilst >>; n; is of course unknown. 

We now turn to the problem of estimating the parameters of the 
distribution which has given rise to the counts x. We give the proof 
that the procedure described in the above example yields the maximum 
likelihood solution. We do this for a general distribution depending 
on k parameters 6, (¢ = 1, 2, --- ; k) which must be estimated. 

Let us consider a discrete variate x capable of attaining integral 
values in the two mutually exclusive and exhaustive ranges denoted 
by ‘7’ and ‘j’ with respective probabilities f(z, 0, --- @,) and f(j, 0, --+ 4). 
Tor these we shall use the short notation f(z, @) and f(j, 6) given by 
(3) and (4) denoting by @ the dependence of f on the multi-parameter 
vector 6, --- 6. We have 


+ = 1 (7’) 


where the summations are respectively over the two mutually exclusive 
ranges 7 and j of the variate x. Further we write 


We now draw a random sample from this population and denote by 
n,; the frequencies observed for the set of integers ‘7’. No frequencies 
are available for the set of integers ‘j’. The maximum likelihood 
equations for the estimation of the 0, are then given by 


8) 
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where the subscript ¢ denotes differentiation with respect to 0,. We 
now introduce auxiliary variables n; given by the proportional allocation 
formulas 


n; = nf(j, 6)/[1 — 7()). (10) 


Substituting >> n,; = n, (8) differentiated with regard to @, and (10) 
in the second term of (9) we obtain 


‘fi, 


so that the system of equations (10) and (11) is identical with the set of 
maximum likelihood equations (9). But equations (11) are formally 
identical with the maximum likelihood equations for a ‘complete 
sample’ of size n + n’ with observed ‘cell frequencies’ n; and n; and 
> n; = n, >> n; = n’. Therefore the iterative procedure described 
above, if it converges, must yield a solution 6, = lim... .6, , 1; = 
lim... -”; of equations (10) and (11) and hence the vector 6 will be a 
solution of the maximum likelihood equations (9). When it is known 
that these latter equations (9) can only have a unique solution the 
iterative process must yield this solution. 

It will be appreciated that the present method is particularly suited 
to situations where the maximum likelihood equations for the complete 
sample are considerably simpler to solve than those for the incomplete 
problem, as is the case for the Poisson. Nevertheless the method will 
still be of great help when the maximum likelihood equations for the 
complete sample must be solved by iterative methods. This is shown 
for the negative binomial in Example 4. 


8. Discrete distributions with grouped frequencies (censoring) 


The classical situation of ‘censoring’ arises when a random sample 
of counts, z, has been drawn from a population and for the counts x 
which exceed the points of censorship 7* (x > 7*) only the total number 
of counts N’ is known whilst the precise values of the counts are un- 
known (have been censored). The counts x with x < 7* are all known. 
(See e.g. the example at the end of this section where N = 240, N’ = 128 
and i* = 2.) 

In this situation it is assumed that in repeated sampling the total 
sample size N is fixed whilst the number of counts N’ in the censored 
section is an observed variable, i.e. is known but not fixed from sample 
to sample. Again, in this section, we consider a more general situation 
of ‘censoring’ which we have termed ‘grouped frequencies’: we assume 
that certain of the sample frequencies have been ‘pooled’ or ‘grouped’ 
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so that only the group totals of these frequencies are known. This is 

illustrated in Example 2 below. It will be clear that ‘censored samples’ 

are covered as the special case when all ‘groups’ consist of single fre- 

quency groups except for one (fairly large) group in which all fre- 

P, quencies for the extreme z-values are pooled. The method of estimation 

1 is described in terms of the following example for which the computa- 
tions are shown in Table 2 below. 


4 Example 2. 


Snedecor [1956, p. 478] gives for 106 litters of 8 pigs the frequencies* 
of litters with 0 or 1 or 2 males (V, = 14), with 3 or 4 or 5 males (NV, = 
73) and with 6 or 7 or 8 males (NV; = 19). Assuming that the number 
of males per litter, x, follows a binomial distribution, to estimate the 
sex ratio of such a distribution of pigs, we require the following notation: 
The variate x is capable of attaining integral values which are arranged 
in G groups, g = 1, 2, --- G. The variate values in each group are 
numbered j = 1, 2, --- , and 


f(j, g; 8) = Pr {a = j-th integer in the g-th group} (12) 
FG; 0) = D0, 1G, 9; 8) = Pr {x in g-th group} (13) 
N, = total observed frequency for g-th group (14) 
: n;, = c-th estimate of j-th frequency in g-th group. (15) 


The procedure of obtaining the maximum likelihood solution is now 
as follows: 
; (a)) From an inspection of the group frequencies N, estimate their 
partitions into individual frequencies on;, . In the example we 
te split N, = 14 into om. = 1, oms2 = 4, oms = 9 and likewise 
N, = 73 into 23, 26, 24 and N; = 19 into 12, 5, 2. This split 
is guided by the knowledge that the distribution is binomial, 
but very rough guesses are quite serviceable here. 
(bo) Using the on;, compute an estimate of the @, from the maximum 
likelihood equations applicable to a complete sample. In the 
example 


1 1 
oN = 106 (22 + 293 + 123) = .516. 


*Actually Snedecor gives the complete set of frequencies which are here grouped as above for 
purposes of illustration. 
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TABLE 2 
ESTIMATION OF BINOMIAL PROPORTION FROM GROUPED DaTA 
o fil 4 34 0030 
1 1 2 4 2.89 .0257 
2 i 3 9 10.77 .0958 
23 21.04 . 2043 
4 2 2 26 28.04 .2723 
24 23.92 2323 
s F i 12 14.13 .1238 
7 3 2 5 4.30 .0377 
N, 14 73 19 14 73 19 1.0000 
DL xn, | 22 293 123 24.43 294.88 119.44 | From table 
i of binomial 
distribution 
19 = 0.516 20 = .51742 


(a,) Using the initial estimate ,@ compute improved values of the 
individual frequencies from ,n;, = N, f(j, g; 10)/F(g; 16). In 


Example 2 the binomial frequencies f(x; 6) = (1 — 


can either be obtained from tables or by direct computation, so 
that we obtain .n,, = .34, im. = 2.89, y3 = 10.77 and the 
remaining ,n,; as shown in Table 2. 

(b,) Using the ,n,; compute an improved estimate 6 from the maxi- 
mum likelihood equations for complete samples. In Example 2 
20 = (24.43 + 294.88 + 119.44)/(8 K 106) = 0.5174. 


A further cycle not shown in the table yields the maximum likelihood 
solution ,;@ equal to 3 decimal accuracy. The proof that the above 
procedure yields the maximum likelihood solution to the problem is 
almost identical with that in Section 2. The maximum likelihood 
equations for the grouped distribution are 


N.F.(g, )/F(g, = 0 (16) 


where the subscript ¢ denotes differentiation with regard to 6,. Splitting 
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up the N, by the proportional allocation 


nip = 9; 9)/F(Q, 4) (17) 
and noting that (13) implies 
Fg, 0) = fi, 9; 0) (18) 


we may rewrite (16) as 
9; 9)/fGj, g; 0) = 0. (19) 


But equation (19) is formally identical with the maximum likelihood 
equations for a ‘complete sample’ of N = )-,, n;, observed frequencies. 


Thus the iteration process dyboa,b, --- , if it converges, will yield a 
solution of (17) and (19) and hence of the maximum likelihood equa- 
tions (16). 


We conclude this section by giving an example of a fairly drastically 
censored distribution to illustrate that the convergence of the procedure 
may still be reasonable when more than half of the distribution is 
‘censored’ (i.e. pooled in one group). 


Example 3. 

Bliss [1953] quoting data of Jones, Mollison, and Quenouille [1948] 
on microscopic counts of soil bacteria gives (Table 4, p. 188) the follow- 
ing distribution of number of colonies per field: 


x = colonies perfield 0 1 2 8 4 5 6+ Total 
n, = frequency ll 37 64 55 37 24 12 240. 


In this example only the extreme tail of the distribution x > 6 is ‘pooled’ 
for the x’-test, but the Poisson distribution was apparently fitted to the 
complete data. We repeat here the fit of the Poisson distribution 
drastically censored at x = 3, i.e. carry out the estimation of the Poisson 
parameter for the distribution 


= colonies per fiedd 0 1 2 3+ Total 
N, = group frequency 11 37 64 128 240. 


The work is set out in Table 3. There are 4 ‘groups,’ the first 3(g = 
1, 2, 3) consisting of the single z values of x = 0, = 1, x = 2 respectively 
and the fourth (g = 4) of the tail x > 3. The initial estimates on,; 
were originally taken as 50, 40, 20, 10, 5, 1, 0; seeing that this total is 
still short by 2 of the required total of N, = 128, the largest frequency 
was increased to 52. This turned out to be a rather lucky ‘guess’ as 
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TABLE 3 
Maximum LIKELIHOOD FIT TO POISSON WITH CENSORED TAIL-SUM 
sg Ng ona; 114; 19) 29) 
01 1 11 
4 37 
33 1 64 
Ss fF 1 52 52.65 52.38 2230 2233 7] 
4 2 40 37.54 37.51 1590 1599 
5 3 20 21.39 21.49 0906 0916 | From 
6 4) 4 10 10.18 10.27 .0431 .0438 | table 
7 5 5 4.16 4.20 .0176 .0179 | of 
8 6 1 1.46 1.50 .0062 .0064 | Poisson 
9 7 0 .47 .49 .0020 .0021 | distribution 
10 8 12 14 .0005 0006 
02 02 .0001 0001 
Totals 112 128 127.99 128.00 .5421 .5457 
DLzn | 165 519 522.59 523.68 
A) 2.850 2.8650 2.8695 


the convergence is established after 2 steps. In Section 5 we shall 
estimate the variance of our estimate as .0155 [see equation (30)]. 
Finally, attention is drawn here to Bliss [1948] dealing with this special 
case of a Poisson distribution with a single censored tail frequency. 

It is well known that the convergence of an iteration can be accel- 
erated by the judicious choice of starting values and in many situations 
special short-cut methods are advocated for their computation. Never- 
‘theless, even if these precautions are taken, convergence may be slow 
in certain problems. We shall therefore describe here a process of 
accelerating the convergence of the present iteration process known to 
computers under the name of ‘approach to the geometric limit.’ This 
method is usually found to be effective when the convergence of the 
iteration process is approximately like that of a geometric series. We 
shall illustrate the use of this procedure for two purposes: 

-@ The conversion of poor starting values to much improved values. 

(ii) The direct computation of the maximum likelihood estimates 

in a poorly converging process. 

(i) Suppose that in Example 1 we had made a rather poor choice of 
starting values as oto = 15, on, = 30, shown in Table 4 below. 
Carrying out the cycles (ao)(be), (a,)(b:) and (a2) as described 
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in Section 1 but carrying fewer figures we reach the following 
three consecutive values ,@ = 2.512, .6@ = 2.761, 36 = 2.900 
with first differences 6,.,,.) = 2.761 — 2.512 = .249, 6.4,2) = 
2.900 — 2.761 = .139. The ratio, g, of the second to the first 
difference is given by 


q = 5063/2) / 8101/2) => .139/.249 = (20) 


and, if we assume that a continuation of the process would 
generate a geometric series with this value of gq, we would 
expect to reach as a limit 


39 + = 2.900 + (.139)(.558) /.442 = 3.075. (21) 


Since the geometric progression of the differences 6 cannot be 
relied upon with certainty it will be necessary to start another . 
cycle of the iteration with a starting value of 6* = 3.075 and to 
carry this until either a limit is reached or until after 2} cycles 
the final value can be computed by a second geometric-series 
projection. In Table 4 the value of ,6 = 3.028 is reached and 
a geometric projection from here yields 3.024. This example 
illustrates that even with poor starting values two sets of 2} 
cycles may yield the limit. 

(ii) If good starting values were chosen but the convergence of the end 
figures is tedious (as in Table 3) one would proceed to the geo- 
metric limit after 23 cycles of the iteration process would have 
have yielded ,@.6 and ;6. In Table 1 we would obtain imme- 
diately 5, (1/2) = —.016, = —.007; q = .44 and 6* = 3.029 — 
.007 (.44)/.56 = 3.023. The method should be used with caution 
for purpose (i) and where convenient short-cut methods are 
available for the computation of initial values, these will be 
usually preferable. The method is, however, rather effective 
for (ii) particularly for improving the end figures of the maximum 
likelihood estimator. (See the 3 computations in Table 5, 
Section 6.) 


5. Computation of variance estimates for the maximum likelihood estimators 


Variance and covariance estimates of maximum likelihood estimators 
6 are usually computed from the expectations of the second derivatives 
E|L,,(0)] of the likelihood function L(@) by substituting the estimates ,6. 
This computation is often facilitated by special aid tables of the elements 
of the information matrix and such tables have in fact been prepared 
for certain isolated cases of maximum likelihood estimation from in- 
complete data. Since situations of ‘incomplete data’ are very varied it 
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does not appear to be practical to provide such aid tables for all such 
situations. It is therefore preferable to use in these situations a numer- 
ical method of estimation based on finite difference calculus and closely 
related to one sometimes used by R. A. Fisher [see e.g. 1953]. This 
appears to be particularly suited to situations in which the maximum 
likelihood equations themselves must be solved by iterative methods 
requiring repeated evaluation of the first derivatives L, of the empirical 
likelihood function (called ‘scores’ by R. A. Fisher). This method can 
be adapted so that it can be used with the present iterative method 
described in the above sections and it will be shown that these iterative 
computations provide, with little extra work, estimates of the variances 
and covariances. 

We first deal with problems depending on a single parameter @ and 
in the next section with the case of two parameters. Starting with the 
case of missing frequencies of Section 2 we note that the first derivative(s) 
of the likelihood function L,(6) with regard to @, are given by the left- 
hand sides of (9) and hence, using auxiliary quantities defined by (10), 
by the left-hand sides of (11). When the maximum likelihood estimates 
@ = 6 are substituted for @ we have of course L,(6) = 0, but when a 
preliminary estimate (say ,@) is substituted we shall obtain a value 
L,(,0) which will be ~0 in all situations in which the maximum likeli- 
hood equation has a unique solution. We can therefore compute an 
estimate of Ly»(6) from the first order divided difference 


Lae = L,(.8)/(.0 6) (22) 


a quantity called ‘rate of change of score’ by Fisher. Higher order 
divided differences may, of course, be used for checking and adjusting 
this computation but, since it is an estimate only, the resulting differ- 
ences are usually small compared with the error involved in the estimate. 

As a check it will usually be wise to compute Lg, from (22) twice 
using two of the trial values .@ and compare the results. Finally, we 
have an estimate of variance in the form 


Var (6) = —1/Ly. (23) 


It should be noted that this estimate of variance introduced by Fisher 
differs from one used on other occasions, viz. 


Var (8) = —1/Lu(6). (23’) 


In the latter formula (23’) the maximum likelihood estimate 6 is sub- 
stituted in a mathematical formula for 1/Z~. , whilst in the former 
formula (23) 6 and a neighboring value .@ are substituted in the formula 
for L; . The second differential Ly, is then estimated by the finite 
difference formula (22). 
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The writer is not aware of any theoretical small sample work on the 
comparison of the relative efficiency of these estimators, either for 
complete or incomplete samples. 

To illustrate the use of (22) we apply it to our examples. For the 
Poisson distribution of Example 1 (see Table 1) we have for the ‘score’ 


L(,0) = (n + yn) — in; + > jin;)/,0 


= 96.499 — (279 + 13.935)/3.052 = .518 


(24) 


and hence 
= (.518)/.026 = 20 
or 
Var 6 = .050. 


A check value computed from ,@ = 2.512 in Table 4 yields Var 6 = .047. 
For an estimate 6 for a complete sample of 96 counts we would have a 
variance estimate of 6/96 = .031. Turning now to the case of pooled 
frequencies in Section 3 it is easy to show by an argument almost 
identical to the one given above that equations (22) and (23) still hold 
for this problem. 

In Example 2 we have for the binomial distribution fitted to grouped 
data 


1 
L,(8) = n,;/0) — nN (25) 
and hence we obtain with ,@ = .51600 and 6 = .51742 
L,(,0) = 751 (887 8 X 106 = 4.82 (26) 


and hence from (22) 
= 4.82/(— 00142) = —339(0) (27) 
so that the estimated variance of 6 is 
Var (6) = —1/Ly0(6) = .000295. (28) 


For a complete binomial sample of N X n = 8 X 106 pigs with sex 
ratio @ = 0.5176 the variance formula is 


(.5176)(.4824)/8 X 106 = .000294. (29) 


Although both the above figures are estimates of variances only, they 
indicate that there is no great loss of information through grouping in 


= 4 
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this example. In the heavily censored Poisson distribution of Example 3 
we find for the trial value, ,@6 = 2.850. 
L,(,0) =n — (do 2n,; + 2N,)/10 
240 — (522.59 + 165)/2.850 = — 1.260, 
Lop = —1.260/.0195 = —64.6, 
Var (6) = —1/Ly, = .0155. 


(30) 


The estimated variance for a complete sample of size 240 from a Poisson 
with estimated mean 6 = 2.87 is 2.87/240 = .0119. 


6. An example of two parameter estimation with truncation* 


As a final example we apply our method to a distribution depending 
on two parameters obtaining their estimates together with estimates 
of their variances and covariances. The example is the negative bi- 
nomial distribution with missing 0-class. 


Example 4. 


Sampford [1955] quoting data of a chromosome breakage study by 
kK. C. Ford gives the distribution of 32 cells with regard to the number 
of breaks per cell. Since the susceptible cells which do not show any 
breaks are not distinguishable from the cells not susceptible to breaks, 
the frequency of the zero-class is not available. We therefore have the 
truncated distribution shown in the third column of Table 5, and as 
given by Sampford. We now use these data to illustrate the fit of the 
negative binomial distribution 


(k— (14+ p) 


*with unknown parameters k and p. Bliss and Fisher [1953] have de- 
scribed a most convenient method of estimating k and p from a complete 
sample and this will form the basis of our procedure. For a complete 
sample of size n + n) = N the first derivatives of the likelihood function 
(called scores by Fisher) are given by 


f(x; k, p) = 


kta (x = 1,2, oes) (31) 


K(k, p) = —* > aps (n + no) In (1 + p) (32) 
and 
Ptk, p) = al. = No_ — kp) (33) 


*This section requires a knowledge of elementary finite difference calculus. 
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where 


A, = nN; (34) 


t=z+1 


is the cumulative sum of the observed frequencies shown in the fourth 
column of Table 5 and ¢ is the mean of the complete distribution of 


TABLE 5 
Maximum LIKELIHOOD FIT TO NEGATIVE BINOMIAL WITH ZERO CLAss MIssING 
Columns (1) (2) (3) (4) (5) (6) (7) (8) 
x § As No ep = flo) =(1+p)* 
0 0 32 40 4.583 . 5637 
1.344 
1 1 41.344 4.499 . 5665 
2 2 6 15 41.818 
4 il 258 
4 4 5 6 | 42.076 4.4549 56808} 
42.088 
5 5 0 6 . 
20 4.231 4373 
f 7 0 5 4.869 
24.869 3.869 4532 
8 8 2 3 1.653 
26. 522 3.759 
.850 
q 27.372 3.7055 46100 
ll 11 1 1 | 27.369 
12 12 0 
13 13 1 0 } 11 2.558 . 2811 
1.512 
12.512 2.471 . 2881 
.438 k=1 
dX in, = 110 12.950 
.179 
13.129 2.4375 29091 
13.128 ] 


a: 
— 
er 
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cm 
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n + mo frequencies. For the truncated problem we must find the roots 
k, p and np, of the equations dL/dk = 0, dL/dp = 0 and 


no = nf(0;k, p)/[1 — k, (35) 


We proceed to solve these equations as follows: three computationally 
convenient ‘pivotal values’ k = 1, k = 3, k = 3 are chosen and for each 
of these values only the two equations dL/dp = 0 and (35) are solved 
by our iterative procedure as shown in Table 5. We describe the method 
for k = 3, note that the equation 0dL/dp = 0 is equivalent to 


p = (36) 
and choose a trial value 92) = 40 as shown in column (5) to compute 
ip = = in,/(n + ond)k = 110/(32 + = 4.583 (37) 


shown in column (7). Next we compute the zero-class frequency for the 
negative binomial from the equation 


k, sp) = 1/0 + = 1/-V5.583 = .5637 (38) 


by using a table of cube roots. This is shown in column (8). Finally 
we obtain an improved value for no as 


imo = nf(O; k, — k, = 41.344 (39) 


which brings us back to column (5). After 23 cycles the geometric 
projection (see Section 4) of om = 40, ym = 41.344, 2no = 41.818 is 
computed from 


= Mo + 4101/2) 9/(1 — = 42.076 (40) 
with g = 81(1/2)/61/2 = .474/1.344 obtained from the first differences of 
the .m) shown in column (6). The projected value ;n) = 42.076 is then 
confirmed by one final cycle of the computations. This process is 
repeated for k = 3 and k = 1 and yields the three maximum likelihood 
estimates p(k) of p for the three pivotal values of k = 3, 3, 1, viz. 


= 4.4549, p(3) = 3.7055, (1) = 2.4375. (41) 


These are set out in Table 6 together with their differences along with 
k and the auxiliary, equidistant tabular argument u = 1/k. Also | 
shown in Table 6 are the three scores 0L/dk = K[k, #(k)] computed 
from (32). We are now ready to find the final maximum likelihood 
solutions k and p. To do this we must find the value of k for which 
dL/dk = 0. We find by inverse interpolation in the table of K that 
@ = 1/k = 2.0287 and hence 


k& = .4929 (42) 
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and by direct interpolation in the table of p 
- p = 3.730. (43) 


These solutions agree with Sampford’s (p. 68) obtained by a method 
requiring a special aid table. Estimates of the variances and covariances 
of & and # can now be obtained immediately from Table 6. Since the 
variances and covariances are computed from the elements of the 
inverse of the information matrix 


_ dk ap (44) 
dk Op ap” 


they can be directly estimated as the partial derivatives of k and p with 
regard to the scores K = 0L/dk and P = dL/dp. We have 


Var = OK? (45) 

Var p= (46) 
ak 
Gy (hp) = = (47) 


A method analogous to equation (45) has, in fact, been used by Bliss 

and Fisher for the variance estimation from the complete sample. 
The partial derivatives 0k/dK and dp/dK can now be estimated immedi- 
ately from Table 6. ‘Since the i(k) satisfy the equation P = aL/dp = 0 
4 the table of the score K as a function of k represents a relation for which 
P is held constant at P = 0. Hence the partial dk/dK can be directly 
computed from the differences of this table as follows: 


(at wu = 1.5) = 1.0327 
a read off from the differences 
ou 
. 
gm (at wu = 2.0287) = sv) by linear interpolation. 
Hence 


(at = 2.0287) = —(h)’ = (.4929)’/.872 = —.278. (49) 


MAXIMUM LIKELIHOOD ESTIMATION 193 
Likewise 


(at u = 2.0287) = (at wu = 2.0287) du (50) 


= .994/.872 = 1.14. 


Finally we obtain the partial dp/dP by using the section for k = 4 
Table 5. The difference between the trial value .p = 3.869 and the 
final value #(}) = 3.7055 can be approximately represented by the 
first order Taylor expansion 


a 
— = + AK. (51) 
Here 
aL 


so that, using (33), we have 


32 + 26.522 


AP = 73 


(3.759 — 3.869) = 


| 


(53) 


Likewise we find 


aL 


or, using (32), we write 


AK = —(32 + 26.522) In 4.869 + (32 + 27.369) In 4.7055 
= — 92.6333 + 91.9466 
= — .687. (55) 


Therefore equation (51) can be written as 
1635 = aP (.171) = aK (.687). (56) 


Equation (56) is a relation between dp/0P and dp/dK at P = 0, K = 
—.0239. Using the same relation* at P = 0, K = 0 we obtain dp/aP 
by substituting the value of dp/@K = 1.14 from (50). We therefore 
obtain 


_ _ 
5.52 (57) 


so that the variance-covariance estimates are given by 
var k= 278, kp = -1.14, Var p= 5.52. (58) 


*In the present example & = .4929 is very close to the pivotal value of k = !9; in general this 
relation would have to be obtained by interpolation between two corresponding relations set up at 
the two pivotal & values neighboring k. 
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Sampford, using the rather complex formulas for the expectation of the 
information matrix [44, p. 65] and inverting this matrix reaches values 
[p. 68] equivalent to 


var k= .276, kp = -1.06, var p= 4.95. (59) 


In comparing (58) and (59) it should be remembered that these are 
values resulting from different methods of estimation. 

Two final remarks are appropriate: 

The pivotal values of k should be simple integers or simple fractions 
so as to facilitate the evaluation of (1 + p)“. They should also be so 
chosen that they cover the maximum likelihood root k. To guide this 
choice some may consider it necessary first to compute a trial k by a 
short-cut method. This may, however, not be necessary as the iterative 
process in Table 5 is fast. One would therefore normally compute 
p(k) for k = 1 and then judge from the sign of 0L/0K (required in 
Table 6) whether one should use further pivotal values of k > 1 or 
k<1. 

The method described in the above example of the negative binomial 
is directly applicable to other two parameter situations. One would 
choose the simpler of the two maximum likelihood equations and solve 
it for the more convenient parameter (0L/dp = 0 is solved for p in the 
above example). This would be done for (say 3) selected pivotal 
values of the other parameter (k in the above example). Such a pro- 
cedure would automatically produce a table of the second score (0L/dK) 
for such values for which the first score (0L/dp) is zero. When this 
procedure is used the simple method of variance-covariance estimation 
is applicable. 
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ON THE PROBABILITY OF SURVIVAL OF BACTERIA IN 
SEA WATER* 


K, Harris 
Public Health Service, Cincinnati, Ohio, U.S.A. 


INTRODUCTION 


The observations analyzed in this paper have been drawn from a 
much larger body of data obtained by personnel of the Public Health 
Service at the Shellfish Sanitation Laboratory, Pensacola, Florida. 
These experiments are concerned with the survival of coliform and 
Salmonella Schottmuelleri organisms added in controlled amounts to 
carboys containing natural sea water at various levels of temperature 
and salinity. 

The primary objective of this work is to determine whether the 
coliform group, nonpathogenic but largely of fecal origin, may be used 
with confidence to indicate the concentration of enteric pathogens such 
as S. Schottmuelleri when both types are being discharged into coastal 
waters. This objective may be satisfied by fitting a general survivorship 
curve to the replicate survival “runs” under each set of environmental 
conditions and analyzing differences in the estimates of the parameters. 

However, within any given set of conditions (e.g., Salmonella in 
water of summer temperatures and high salinity) individual runs 
showed relatively wide variation in estimated bacterial densities, 
particularly at times of one or more days after the start of an experiment. 
Realizing, then, that average survival curves, although useful for 
comparison purposes, would be limited in their ability to predict actual 
densities, we were led to formulate a second objective. Assuming that 
knowledge of the paths and velocities of off-shore currents allows us to 
predict the time required to traverse the distance from some point of 
discharge, e.g., a sewage outfall, to some other location such as a potential 
swimming or shellfish growing area, we would like to be able to estimate 
the probability that X or more organisms out of a given initial concen- 
tration will survive this time interval. 


*Presented at a joint session of the Biometric Society and the Institute of Mathematical Sta- 
tistics, annual meeting of the American Statistical A iation, September 10, 1957, Atlantic City, 
New Jersey. 
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*Geometric mean of the two middle-sized MPNs. 
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The present paper is devoted entirely to this second objective and 
proceeds by constructing and fitting a suitable probability model. 
Of course, in this area of environmental sanitation no model based on 
laboratory data can accurately predict events under typical natural 
conditions. Nevertheless, for reasons given in the final section below, 
results obtained from the theory developed here may prove useful 
without further modification. 


DATA AND ANALYSIS 


The procedure most commonly followed to estimate bacterial 
densities in water is the dilution test involving, here, three decimal 
dilutions using five tubes per dilution. From the number of sterile 
(or fertile) tubes in each set of five, we are led under certain assumptions 
to a maximum likelihood estimate of bacterial density per ml., the 
so-called MPN or “Most Probable Number.”’ The underlying theory 
has been summarized by Cochran [1950]. 

To show more clearly the wide variation observed in MPNs at 
given sampling times, and to provide illustrative material for analysis, 
Table 1 has been prepared, listing the complete data—including numbers 
of sterile tubes at each dilution—from replicate runs with Salmonella 
in sea water of relatively high salinity (20-30 ppm) at 25-29°C. 


Compound Binomial Model 


Let 6, denote the density of organisms per ml. at time ¢ (in days) 
after the start of any survival run. Then, variation in replicate MPNs 
is assumed to arise from inherent variability in the MPN estimate of 
a given density 6, , compounded with variation in 6, according to some 
probability model. Variability in 6, may stem from the random effects 
of many uncontrollable factors, e.g., varying plankton densities, or 
different concentrations of salts and nutrients in the sea water medium. 
Variation in initial loadings, 5, , would also be reflected in 6, although 
the effect of this source of variation would decline rapidly as ¢ increased. 

The model developed below applies to the case of an MPN estimate 
based on a single dilution. Although the mathematics is not difficult 
in this case, the writer has not been able to generalize results to the 
multidilution case. In these experiments, however, much useful and 
reliable information may be gained by studying the data from a single, 
carefully selected dilution at each sampling time. 

Let f(s) be the probability of s sterile tubes when n tubes are planted, 
each with v ml. of liquid containing an average of 6 organisms per ml. 
Then, we suppose that 

jo) = (Rea aa, 


t 
= 
} 
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where 
[1950)]). 

Thomas [1952] has suggested the Gamma distribution as a reasonable 
form for ¢(6) and has presented evidence, based on plate count data, 
to substantiate this model. Following his lead, we write 


= (") if — di 


is the probability of a sterile tube, given 6 (cf. Cochran 


pti 1 
= (n) f (1 — 6°)" (log 1/6)” dé, 
p! Jo 
where @ = 
Expanding the binomial (1 — 6")"”* and integrating term-by-term, 


n 


pti 1 1 
n-s8 1 
+( 2 (1) 


n n-s q 
From this formula, a general equation for factorial moments is 
easily derived. In fact, any model of ¢(6) which permits expansion of 
the binomial (1 — e~*’)"”* to be followed by term-by-term integration 
or summation (e.g., the normal, Poisson, Gamma or any simpler form) 
would lead to an expression for f(s) of the general form, 


whose j-th factorial moment is given by the equation, 


[See appendix for derivation of Eq. (2)]. Applying Eq. (2) to (1), 


Ke) = nf: + -) n(; + 7 


and 


(3) 


where m = q/v. 


| 
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To use this theory in multidilution experiments, one dilution must 
be selected which, if it had been the only one tested, would have yielded 
an MPN estimate closest to that actually produced by the multiple 
dilutions. In a three-dilution test where the volumes are in decimal 
ratio, this selection is usually not difficult since one of the dilutions 
bears the main burden of estimation, greater dilutions yielding all or 
almost all sterile tubes, and smaller dilutions practically none. This 
is brought out in Table 1, where the critical volumes, day by day, are 
fairly obvious: .1, .1, .1, 1, 1, 10 and 10 ml., respectively. 

Applying Eqs. (3) to the distribution of sterile tubes at each of 
these volumes (the subscript ¢ for time is understood in the following), let 


mean 8; 


S, = mean s(s— 1). (A) 


n(n 
Then, 


K = log S,/log S, is an estimate of log / 


Thomas [1952] has tabled 2 (called a in his table) as a function of K. 
Finally, 


wen _ bt! 
Ka) = (5) 
Estimated 
Var 6 (om)? 


K lies in the range 1 to 2, approaching the lower limit as m — 0 
and the upper limit asm — . An observed value of K greater than 2 
implies that Var 6 = 0, $(6) collapsing to a single point; that is, the 

variation in s is no more than would be expected in repeated sampling 
from a constant density, 6. In such a case, 6 is estimated from the 
usual MPN formula: 


Table 2 lists results obtained on fitting Eq. (1) to the distributions 
of sterile tubes at selected volumes within the experiments presented 


1 
§ = ogi 
v 
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in Table 1. Agreement with the model is generally good at each sampling 
time as well as when total frequencies over all sampling times are 
examined. The model exhibits considerable versatility, being able to 
produce bimodal, U-shaped curves and uniform distributions (m = 
p + 1 = 1) as well as typical humped forms. 


TABLE 2 


OBSERVED FREQUENCIES OF STERILE TUBES COMPARED WITH EXPECTED 
FREQUENCIES UNDER GAMMA DisTRIBUTION MODEL (gq. 1), aT Eacu 
SAMPLING TIME AND OvER ALL SAMPLING TIMES 


Sampling Time 


0 days 3 day 1 day 14 days 
(.1 ml.) (.1 ml.) (.1 ml.) (1 ml.) 


s Obs. Exp. Obs. Exp. Obs. Exp. Obs. Exp. 


0 1 .87 1 .87 0 .25 5 5.02 
1 2 2.08 2 2.08 1 56 2 Bee 
2 2 2.37 2 2.37 1 94 0 .76 
3 2 1.73 2 1.73 1 1.37 1 .34 
4 1 .78 1 .78 2 1.96 0 .14 
5 0 17 0 Bf 3 2.92 0 .04 

2 days 3 days 4 days 

(1 ml.) (10 ml.) (10 ml.) Total 


0 2 1.23 1 1.30 1 74 11 10.28 
1 0 1.26 2 1.60 0 .63 9 9.92 
2 1 1.26 2 1.65 1 64 9 9.99 
3 2 1.32 1 1.48 1 .82 10 8.79 
4 2 1.38 1 1.21 1 1.19 8 7.44 
5 1 1.55 1 76 4 3.98 9 9.59 


An index to the goodness of fit is provided by the likelihood ratio test 
rz r I] 


where r denotes the total number of replicate experiments, r,, the 
number yielding s sterile tubes at sampling time ¢, and f,(s) the prob- 
ability of s under the model fitted at ¢. 


| 
: 
| | | 
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ee a s | Obs. Exp. Obs. Exp. | Obs. Exp. |Obs. Exp. 
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The distribution of —2 log, \ is closely approximated by x” even 
for small values of r. Degrees of freedom here equal n — 2, or 3. Since 
the observed distributions of sterile tubes at the different sampling 
times are mutually independent the test criterion has been x37 = 
—2 557 log, \+ , where T is the number of sampling times during a 
given set of experiments. Calculations for each of eight sets involved 
in these studies have yielded values of x’/d.f. ranging from .48 (for 
the set presented here) to 1.42, five of the values lying between .9 and 
1.1, with an over-all average of .96. 

Estimates of E(é,) and Var 6, from Eqs. (6) are shown in Table 3. 
At t > 3 day, the coefficient of variation of 6, was generally large but 
independent of £(6,). When t < } day the organisms had not yet 
become affected by the numerous sources of variation in the water and 
Var 6, was relatively small. 


TABLE 3 
Estimates or MEAN DENsity PER ML. [/(6,)] AND 


VARIANCE OF DeEnsity (Var 


Time 
(in days) Var 
0 9.68 11.00 
3 9.68 11.00 
1 3.93 18.72 
13 3.03 3.12 
2 .943 .993 
3 .102 .00624 
4 .059 01204 


Variances and efficiency of moment estimates 


Approximate variances and covariance of the moment estimates 
m and p are given by the equations: 


Var m D Var S, D 


( 
Cov (S,, 8.) + ap + 1 Var S, 


a(p + 1) 
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J” Ver p = ( Var 8S, Cov (S, , 


(6) 


J* Cov (mn, p) am ‘ap + 1) Var S, + am ‘ap + 1) 


4 Cov (8, , & Var S., 


am om Ap + 1) 
where 
dm Ap+1) dm Apt 1) 
1 m 
1 m 
= n(n 1) E|s(s 1)] (; + ° 
Under the model, 
Var S, 
m p+1 m p+1 p+1 
= - +() [1 2 } 
Var 


Cov (S, S,) 


1 m m \'! m 


where r denotes the number of replicate experiments. The estimates 
m and # are used in place of the true values. 

The efficiency of these moment estimates relative to maximum 
likelihood estimates has been examined. Here, 


é 

1 

+(—“) Var S, 
om = 

| 
4 
n 
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and, letting 


_ _ — 8 m 


the m.1. equations are 


om m 
dlogL _ 


apt)” 
which can be solved ouly by a very laborious process of successive 
iteration. 

The determinant of the asymptotic m.l. information matrix is 

(p+ Ap ( 

while the determinant of the variance-covariance matrix of the moment 
estimates reduces to [Var S,-Var S, — (Cov S, , S2)*]/J’. 

The product of these two expressions is the reciprocal of the efficiency 
of the moment estimates. 

Efficiencies have been calculated for the following forms of the 
density function, f(s): U-shaped (m = .01, p + 1 = .15), bimodal but 
asymmetric (m = p + 1 = .3), uniform (m = p + 1 = 1), reverse 
J-shaped (m = 2.1, p + 1 = 4.5), humped (m = p+ 1 = 8), and 
J-shaped (m = 5.4, p + 1 = 1.13). In each of these cases, representing 
a wide range of values of the parameters, the moment method of estima- 
tion was found to be at least 95 per cent as efficient as maximum likeli- 
hood. 

A recent paper by Sampford [1955] was helpful in preparing this 
section on variances and efficiency. 


APPLICATIONS 


Application of this probability model rests on the following unrealistic 
assumption, namely, that a dose of bacteria discharged into a natural 
body of water moves with the currents as a packet whose density 
depends solely on the organisms’ survival pattern. This assumption 
is necessary since the present experiments have been conducted under 
stationary conditions in laboratory carboys. In natural waters longi- 
tudinal mixing (at least) occurs so that the time spent by an organism 
in traveling a certain distance is itself a random variable whose distri- 


dlogL r 
8 
\ 
pad 
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bution must be combined with the distribution of the proportion of 
organisms surviving to that time before an approximate model of the 
real situation may be attained. 

In addition, the dilution provided by a large volume of natural, 
uncontaminated sea water has not been considered. Dilution varies 
with the tides and turbulence of the waters although an average dilution 
factor would not be difficult to obtain for any given area. 

Mixing and dilution operate to reduce the density of bacteria 
arriving at some location to an amount less than would be predicted 
from these laboratory experiments. The model in its present form 
yields estimates of bacterial survival which might occur in nature only 
under the most unfavorable conditions. These estimates therefore 
indicate maximum hazard and as such may prove useful in planning 
and administering programs of sanitary control. 

We wish to estimate the probability that, given a large initial dose 
of organisms, the concentration surviving the time necessary to traverse 
some critical distance (e.g., from sewage outfall to potential shellfish 
growing area) will exceed some much smaller number. The appropriate 
variable is 6,/5, . To simplify matters, we select a relatively accurate 
estimate of 5, —say, the median of replicate multidilution MPNs— 
and assume that 6,/6, is distributed approximately in a Gamma form 
with parameters p, and g,59 = v,m,6,. For the Salmonella data given 
in Table 1, 4, is set equal to 10.1 organisms per ml., the median MPN 
at zero time. 

Consider, for example, survival probabilities after three days’ 
traveling time of a heavy outfall of Salmonella organisms, say, 10,000 
per ml. The relevant statistics obtained in these experiments under 
conditions of high salinity and temperature are 


Ds = .68, V3Mls = 10 x 1.64, 
leading to the estimate 
4 
= X 10, 
or 101.4 organisms per ml. 
Using Eqs. (6), Var #; = 2.3251, Var m3 = 3.3241, and Cov (3 , ps) 
= 2.6176. The approximate standard error of E(é;) is given by the 
expression 


A 2 
10" Var ps + Var tis 


3 
= o(& +1) Cov (m;, | 
3 
= -+40.5 organisms per ml. 
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Pearson’s Tables of the Incomplete Gamma Function [1934] are 
entered with the arguments f; = .68 and u; = (vst36o/ Vp, + 1) y 
= 127.8ly, seeking the probability that a proportion y or greater of 
the initial concentration will survive. Thus, setting y = .02, or 6; = 200 
organisms per ml. from an initial loading of 10,000 per ml., we estimate 
from the tables that 11 per cent of the time the actual concentration at 
this distance (three days’ traveling time) would exceed 200 per ml. 
About 70 per cent of the time the density would exceed 50 organisms 
per ml. 

Upper confidence limits for these percentages, or an upper tolerance 
limit for the concentration of organisms exceeded in a given percentage 
of cases, may be approximated by calculating a joint confidence region 
for m and p and selecting a pair of values which minimize u. This is 
not easily done when m and p have been based on small numbers of 
observations, as here, since these estimates will not be bivariate normal 
and hence their quadratic form cannot be used to obtain a joint con- 
fidence region. We may rely, however, on the x” approximation to 
the likelihood ratio test and attempt to minimize u’ = m/Vp+1 
subject to the restriction that —2 log, \ equal, say, 7.81, the 5 per cent 
level of x3y.;. . In the present case, an analytical solution does not 
seem possible; however, by moving from*the estimates in the direction 
of smaller m and larger p, we can arrive at a trial-and-error solution. 
For the example here, m; = 1 and p; = 1.40 yield —2 log, \ = 7.79 
and uz; = 65.19y. Entering Pearson’s tables with these values of u 
and p, we conclude with approximately 95 per cent confidence that the 
3-day concentration will exceed 200 organisms per ml. no more than 
52 per cent of the time and 500 per ml. no more than 7 per cent of the | 
time. 

Suppose we desire the concentration at this distance to be no greater 
than one organism per ml. more than 5 per cent of the time. Using 
ms and pj; we estimate from Pearson’s tables that the initial concentra- 
tion should not exceed 40 organisms per ml. A lower tolerance limit for 
this density, based on m; = 1 and p, = 1.40 is found to be 19 organisms 
per ml. 

These are, of course, purely illustrative examples to show the possible 
application of this probability method to real problems in the sanitary 
control of coastal waters. Analysis of the mean survivorship curves 
mentioned in the introduction is likely to find that variation in certain 
factors has not produced significant change in the over-all survival 
pattern. In this case some of our results may be averaged and more 
precise estimates obtained. 
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APPENDIX 


Derivation of Eq. (2): 
Starting from the general form, 


= (") Jae + 0, 
the j-th factorial moment of s is given by 


k=0 


s!kiin —s — k)! 

Rearrange the terms in this summation into groups of sums such 
that, for the terms within a given group, s + k is a fixed quantity, say 
r (not to be confused with r in the text), where r > j since according 
to Eq. (1’) we need not consider any value of s < j. Then, 


Efs(s — 1) --- (8 + 


r! 


| Sat r)! (—1) k\(r g(s +k= 


When r = j, the sole term in brackets is that for k = 0 and s = j, 
namely, [n!/(n — j)!]g(j). Since Eq. (2) in the text states that this 
term itself is equal to the j-th factorial moment, it must be shown that 
for every r,n > r > j, the summation in brackets is zero. Factoring 
out g(s + k = r), we require that >-iz, (—1)*r!/[k\(r — k — 9)], or 
rather, >-7-} (—1)*(r — j)!/[k\(r — k — 9)!], equal zero. But this is an 
alternating sum of binomial coefficients and is well known to equal zero. 
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EFFICIENT ESTIMATION OF THE RELATIONSHIP 
BETWEEN PLOT SIZE AND THE VARIABILITY OF 
CROP YIELDS* 


W. H. Harneway anp E. J. 


Rockefeller Foundation, Agricultural Field Staff, and the Institute of Statistics 
Raleigh, North Carolina, U.S.A. 


1. Introduction 


The optimum size of plot in field experimentation depends on the 
relationship between fixed costs and costs varying with number of units, 
and on soil variability. Perhaps the most useful measure of soil hetero- 
geneity yet devised is that of Smith [1938], who showed empirically 
that the logarithm of the variance between plots of a given size was 
linearly related to the logarithm of the size of the plot. In the present 
paper we consider only the relationship between size and variability. 
The objects of the paper are, firstly, to show how efficient estimates of 
the constants in this relationship may be determined, and secondly, to 
illustrate a general method of determining efficient linear estimates 
when the data are, as in the present instance, correlated and of unequal 
variability. 

Koch and Rigney [1951] demonstrated that the regression coefficient 
of the logarithm of variance on the logarithm of plot size could be 
estimated from experimental data in which treatment effects are present, 
as well as from the data of uniformity trials. They noted that Smith 
had recommended that, in estimating the regression coefficient 8, the 
variances of the different sized plots should be weighted by their respec- 
tive degrees of freedom. In fact, since the variance estimates for different 
size of plot, both in uniformity trials and experimental data, are built 
up from common components, they are frequently highly correlated, 
so that a simple weighting by degrees of freedom is not accurate. Koch 
and Rigney point out this difficulty for experimental data, but do not 
seem to have realized that their arguments apply with equal force to 
uniformity trial data. 


The present paper presents a method of weighting observed variances 
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of different-sized plots which leads to an unbiased estimate of 8 with 
asymptotically minimum variance. It is applicable both to uniformity 
trial data and to experimental data; in the latter case the analysis of 
variance is in effect reconstructed to simulate one derived directly from 
uniformity trial data, in the manner suggested by Koch and Rigney. 


2. Regression Analysis with Correlated Unequal Errors 


In the familiar applications of regression analysis, it is often assumed 
that the values of the dependent variable are uncorrelated and of equal 
variance and that the values of the independent variable are errorless. 
Alternatively, if the values of the independent variable are subject to 
error, it is assumed that the deviations of the dependent variable from 
the hypothetical regression function are uncorrelated and of equal 
variance. The first case is in fact a particular case of the second. In 
either case, the application of the method of least squares to the estima- 
tion of the regression coefficient is straightforward, the quantity to be 
minimized being the sum of squares of the deviations. Thus, if « and y 
are the independent and dependent variables respectively, then the 
sum of squares to be minimized is 


Bx,)° 
leading to the familiar estimates 


b yi(z; — (x; — #)° 


i 


ll 


bi 
for 8 and a respectively. 

A case of less frequent occurrence is that in which the values of y 
are correlated and of unequal variance. Although the method of least 
squares will still give an unbiased estimate of the regression function, 
this method will in general be inefficient, so that the estimates b and a 
will not be the most accurate obtainable from the data. If the variances 
and covariances of the y-values are known, efficient estimates of the 
constants may be determined; and even if the variances and covariances 
are only estimated from data, their use enables estimates with approxi- 
mately minimum variance to be determined. 

In this general case, suppose that there are n y-values, that the 
variance of y; is o;; , and that the covariance of y; and y, is o;,. Then 
it is known (see, for example, Rao [1952] Ch. 3) that, subject to the 
usual assumptions of normality, the quadratic form 


— a — — — Bay) 
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is distributed as x” with n degrees of freedom. Here the terms o™ are 
the elements of the inverse of the covariance matrix of the y; . 

When this quadratic form (which takes the place of the sum of 
squares in the usual analysis) is minimized with respect to @ and 8, 
it yields efficient estimates of the constants as follows: 


and 
a=y%-— bi. 


Here € and g are weighted means; thus 
== =. o'x,/ >, = 
i k i k 


From consideration of b as a linear function in the y; it is readily 
verified that its variance is in fact simply 


The quadratic form now has the distribution of x’? with n — 2 
degrees of freedom, since it has been minimized with respect to two 
parameters. Thus, besides estimating the constants, we may also test 
the significance of the residual form, which indicates departure from 
linear regression. 

When variances and covariances are estimated, the results given 
here still hold approximately. In the following discussion we shall use 
these general results, assuming that the errors in estimation of the 
variances and covariances can safely be ignored. 


3. Estimation from Uniformity Trial Data 


Koch and Rigney showed that a uniformity trial subdivided to 
simulate a split-plot or lattice design could be analyzed in the manner 
shown below; a randomized block arrangement could similarly be 


Degrees of Mean Expectation of 
Source Freedom Square Mean Square 
Replications d—1 Vi S +aP + abQ + abcR 
Blocks within 
replications — 1) V2 S+aP + abQ 
Plots within 
blocks cd(b — 1) V3 S+aP 


Subplots within 
plots bed(a — 1) Ve S 
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superimposed on the trial, though it would not provide so much informa- 
tion about the relationship of variability to plot size. 

We shall denote the variances of plots of various sizes, reduced to a 
subplot basis, by V’. Thus, if a complete replication is regarded as a 
large plot, its variance V{ is equal to the replication mean square as it 
appears in the analysis of variance; 


Now, considering blocks as the next smaller size of plot, we see that 
the variance between blocks contains, in addition to the variation due 
to blocks within replications, that removed by the stratification of 
groups of blocks into replications in the analysis of variance. Thus the 
total sum of squares for blocks is 


and, since there are cd blocks, the mean square is 
Vz = [de — + — — 1). 
Similarly, the variance between plots over the entire area is 
Vs = [ced(b — 1)V3 + de — 1I)V. +  — 1)Vi)/(bed — 1) 

and the variance between subplots over the entire area is 

Vi = [bed(a — 1)V, + cd(b — 1)V, + de — 1)V, 

+ (d — 1)V,]/(abed — 1). 


These formulas are formally identical to those given by Koch and 
Rigney, who expressed their results in terms of components of variance. 
Smith’s regression coefficient 8 is defined by the formula 


E(log V.) = E(log Vi) — B log x 


where x is the number of units per plot, V, is the variance among plots 
one unit in size, and V, is the variance of mean yield per unit area for 
plots of size x units. This formula requires a “finite population” 
correction when the size of block is small compared with the size of plot. 
This correction is discussed in Smith’s original paper. We propose to 
ignore it because the variance estimates which it affects most are those 
least accurately estimated. For purposes of estimating optimum plot 
size, the coefficient 6 is alone of interest. In the computations suggested 
by Koch and Rigney, the values of V, are obtained by dividing each 
value of V’ by the number of units per replication, block, plot or subplot, 
thus putting them on a unit basis. According to them, 8 is estimated 
by means of an unweighted regression coefficient: 


if} 
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yi(xj; — 


where y = log (V’/zx), 2’ = log a, and # = >>; x//n. 


As was pointed out to one of the authors by D. D. Mason, Depart- 
ment of Experimental Statistics, N. C. State College, application of this 
formula for b, sometimes gave results greater than 1, which are un- 
acceptable on physical grounds. It was realized that the above estimate 
(1) would often be inaccurate owing to the equal weighting of y-values 
of differing variability. It was therefore decided to apply to the different 
terms in the sums of squares and products defining the regression 
coefficient, weights that would lead to an estimate with minimum 
variance. The appropriate weights are the elements of the inverse of 
the covariance matrix (i.e. the information matrix) of the values of y. 
If these elements are designated w;, , the estimate is 


U 


where 

The weights w;, will have to be estimated from the data and will 
be to that extent inaccurate; but apart from this source of error, the 
effect of which we do not consider, the estimate will be of minimum 
variance; this variance is in fact 


1 1 


When, as is often the case, there are more than two variance estimates 
from which to compute the regression, we may also test the significance 
of departure from regression. The weighted total sum of squares of 


the Y; is 
where 


with n — 1 degrees of freedom, n being the number of variance estimates. 
The sum of squares attributable to regression on 2} is 


U?/T. 
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Hence the sum of squares for departure from regression is 

V — U’/T, 
which is distributed approximately as x” with n — 2 degrees of freedom, 
and may be tested accordingly. 

It now remains to estimate the weights. Since V,; , V2, V3 , and V, 
are independent, and their variances are 2V{/(d — 1), 2V2/d(e — 1), 
2V;/cd(b — 1), and 2V%/bed(a — 1) respectively, it is not difficult to 
determine the variances and covariances of V{ , Vi , Vj, and Vi, 
which are linear functions of the former set. In fact, not only the 
variance of V{ , but also its covariances with the other V‘ , are propor- 
tional to 

Likewise the variance of V{ is estimated as 


2[d(e — 1)V2 + (d — 1)Vi)/(ed — 


and its covariances with V{ and V{ are proportional to this. Thus we 
find the covariance matrix of the V/ to be as follows: 


D D D D : 
(d—1)° (d—1)ed—1) (d—1)(bed—1) (d—1)(abed—1) 
D C+D C+D C+D 
(d—I(ed—1) (ed—1)(abed—1) 
D C+D B+C+D B+C+D 
(d—1)(bed—1) (ed—1)(bed—1) (bed—1)° (bed —1)(abed— 1) 
D C+D B+C+D A+B+C+D 
| 1)(abed—1) (ed—1)(abed—1) (bed—1)(abed—1) (abed—1)? 


where 
A = 2bed(a — 1)Vi, B = 2cd(b — 1)V3, 
C = 2d(e—1)V;, and D=2d—1)Vi. 


The inverse matrix is found to be even simpler in form; as may be 
verified, it is 


—(d=1)(ed—1) (1,1) —(d—1)(bed—1) 

0 = bed (bea — (ted =D) 
0 0 — (bed —1)(abed— 1) (abed— 1)? 

a A A a 


| | 
of 1 1) —(d—1)(ced—1) 
. (d—1) ( ot D C 0 0 
| 


PLOT SIZE AND CROP YIELDS 213 


The weights for y; (= log V/) are obtained by multiplying each 
row and each column of this inverse matrix by the corresponding V/ . 
This result follows from the approximate formula 


Cov (log Vj , log Vi) ~ Cov (Vj, Vi)/ViVE , 


which results from a routine application of the method of ‘propagation 
of error’. If, as is usual in practical computation, logarithms to base 
10 rather than natural logarithms are taken, the weights will need to be 
multiplied by the factor 


M~ = (logo e)” 
= 5.302. 


We shall deal here with the transformation to natural logarithms, and 
indicate the adjustments necessary for common logarithms below. 

Thus, from the inverse matrix, 

Witty; = ~ , Vi). 
The weights may thus be determined from the inverse matrix without 
too much difficulty. 

It will be found that the sum of the elements of the weight matrix 
is equal to half the total number of degrees of freedom for the sums of 
squares from which variance estimates are derived. This may be seen 
in the following way. If the variances are unaffected by size of plot, 
then all the available sums of squares are estimates of the same basic 
variance. The different estimates of the logarithm of the variance, 
derived from different lines of the analysis of variance, are independent, 
and have asymptotic variance equal to twice the reciprocal of the 
corresponding degrees of freedom. Consequently the information 
from each is half the degrees of freedom, whence the total information 
is half the total degrees of freedom. 

Thus, for data from uniformity trials, the sum of the weights is 


4(abed — 1), 
while for data from split-plot experiments, the sum of the weights will be 
$be(ad — 1). 


Similar results may be derived for lattices and other types of experi- 
mental design. They provide a convenient check on the computation 
of the weights. 

To determine the regression coefficient and to test the departure 
from regression, the calculation is best carried out in stages, as follows. 
Let 


xX, = W 
i 
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and 
Y, = - 

Then the sum of squares of 2’ is 

Similarly the sum of products of y with z’ is 

U= Xiyi — Y)/ wis 

and the sum of squares of y is 


Then b, = —U/T. 


The variance of the estimate b, is, to the degree of approximation of 
the analysis, 


Hence, approximate confidence limits for the population regression 
coefficient 8 are 


b, + 


t being the normal deviate at the required level of probability. 
Departure from linear regression is tested, as indicated above, by 
means of 


V U/T 


which is regarded as x” with n — 2 degrees of freedom. 
When common logarithms are used, the value of b. is determined as 
above, but its variance is now 


T~*'/5.302 = 0.18867", 
and the corresponding confidence limits are 


b, + 0.4343¢T"?. 


The sum of squares for departure from regression is also altered, to 


5.302(V — U?/T). 
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It should be observed that the key to these computations is the 
covariance matrix of the variances V‘ of the plots of different sizes. 
Because these variances are expressed as linear combinations of the 
original mean squares, which are independent, and not in terms of the 
variance components, which are correlated with one another, the result- 
ing covariance matrix, and its inverse, take on a relatively simple form. 


4. Estimation from experimental data 


When variance components are to be estimated from experimental 
data, the estimates are calculated in the same way as from uniformity 
trial data. However, since a number of contrasts are given over to the 
estimation of treatment effects, the different plot and block variances 
are estimated with fewer degrees of freedom, and hence less precision, 
than they could have been in a uniformity trial. Apart from this 
complication, for which allowance must be made in determining the 
weights for the various components, the determination of a linear 
unbiased estimate with asymptotic minimum variance follows the same 
lines as that given in the previous section. The method is illustrated 
by the analysis for a split-plot experiment in the form given by Koch 
and Rigney. It will be noted that, in this model, it is assumed that 
block-treatment interactions do not exist. 


Degrees of Mean Expectation of 
freedom Square mean square 
Replications d-1 Vi S + aP + abQ + abcR 
Treatments (1) c-1 S+aP + abQ+ 
treatment effects 
Error (1) (c — 1)(d — 1) V2 S +aP + abQ 
Total between 
whole plots cd —1 
Treatments (2) and 
interactions c(b — 1) S + aP + treatment effects 
Error (2) — 1) (d — 1) Vs S +aP 
Split-plots cd(b — 1) 
Sampling error bed(a — 1) Vi S 


As for a uniformity trial, the estimated variance of plots the size of 
a complete replication is V, . . Since it is estimated with the full d — 1 
degrees of freedom, its variance is as given in the previous section. 

In estimating the mean square for blocks (i.e. whole plots) we must 
allow for the fact that, of the d(c — 1) contrasts between blocks within 
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replications, only (ce — 1)(d — 1) are available for estimating the variance, 
the other c — 1 containing treatment effects. Thus, as before, the 
estimated variance between blocks is 


= [de — 1) V2 + (d — — 0), 
but its estimated variance is now increased to 
2[d’*(c — 1)V2/(d — 1) + (d — 1)Vi)/(ed — 1)’. 


The variance of V, has similarly to be adjusted by a factor d/(d — 1). 
The analysis now proceeds as for uniformity trials, and the inverse 
matrix is-as given above, provided we redefine 


A = 2bed(a—1)Vi, B= — 1)Vi/(d — 1), 
C = — 1)Vi/(d — 1), and D=2d—-1)Vi. 
5. Numerical example 


The computations required in the proposed method are illustrated 
in a numerical example, the data for which, set out in Table 1, were 
kindly furnished by D. D. Mason. 


TABLE 1 


SoyBEAN YIELD TRIAL ConpuctTep By C. A. Brm, 
U.S. DEPARTMENT OF AGRICULTURE, AT WILLARD, NortTH CARo.ina, 1956. 


Degrees of 
Source Freedom Mean Square 
Replications 2=d-1 452 = DV, 
Varieties ll=c-1 30401 
Experimental Error 22 = (d — 1) (ce — 1) 10589 = V2 
Rows in Plots 36 = cd(b — 1) 5938 = V3; 
Subplots in Rows 72 = bed(a — 1) 2862 = V, 


Here a = 2,6 = 2,c = 12,d = 3. 


In the determination of the weights we can work with multiples of 
the V‘ more conveniently than with the V/ themselves. This device 
makes the computations simpler as well as more accurate. We have 


2V, = 904, 
(ed — 1)V2 = 35V3 = 2V, + 33V, 350341, 
(bed — 1)V3 = 71V3 = 2V, + 33V, + 36V; = 564109, 
(abed — 1)Vi = 148V{ = 2V, + 33V, + 36V; + 72V, = 770178. 
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This gives 

Vi = 452, V3 = 10010, Vi = 7945, and V4 = 5386. 

The number of units per plot corresponding to the different-sized 
plots the variances of which are V{ , V3, Vi, and V{ are 48, 4, 2, and 1 


respectively. Putting the variances V’ on a unit basis and taking 
logarithms, we obtain the values given in Table 2. 


TABLE 2 


LoGARITHMS OF RELATIVE PLoT S1zEs (2’) AND 
Unit VarRIANCEs (7) 


x y 
0.0000 3.7313 
0.3010 3.5891 
0.6021 3.3984 
1.6812 0.9739 


The unweighted regression coefficient is then 


Dery - /| - 


— 1.6697 
+1.7334. 


As Koch and Rigney point out, 8 is an index of soil variability; it 
should vary between zero and unity. A value of zero indicates perfect 
correlation (extreme uniformity) among the units making up a plot; 
a value of unity indicates no correlation. Clearly the estimate b, 
obtained in the present example can have no unambiguous physical 
interpretation. Here it is apparent that weighting a low mean square 
for replications based on only two degrees of freedom, equally with 
others based on many more degrees of freedom, has led to an unreason- 
able estimate of soil variability. 

Using the method proposed in the present paper, y and 2’ are as 
before. The weights are the elements w;, of the information matrix 
of the y. To obtain these numbers it is convenient first to calculate 


b, 


Il 


A = 2bed(a — 1)Vi = 1179.5 xX 10° 
B = 2cd*(b — 1)V3/(d — 1) = 3808.1 X 10° 
= 2d°(c — 1)Vi/(d — 1) = 11100.6 X 10° 

= 2d — 1)Vi = 0.8 x 10° 
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then 
1 


nvir(t + 1) = 1.00 


= Wa = —[(d — — 1)V2]/C = —0.03 


Wis = Wy, = Wis = Wa, = Oz 


The remaining elements are computed in similar fashion. The completed 
information matrix is 


1.00 —0.03 0.00 0.00 ] 
70038 43.29 —51.90 0.00 
0.00 —51.90 353.36 —368.34 
0.00 0.00 —368.34  502.89_ 


the sum of whose elements is }bc(ad — 1) = 60, as may be verified. 
We now compute the set 


Thus 
Y, = (1.00)(0.9739) — (0.03)(3.3984) = 0.87; 
similarly 


Y, = —39.19, Y; = —282.52, Y, = 554.42. >) Y; = 233.58. 


In the same way we compute the X, : 


X, = 1.66, X2.= 10.39, X;3 = 75.11 


X,= = —23.71. 
Then 
T= X — X)?/ Wik 
= 31.65 — (—23.71)’/60 = 22.28. 
Similarly 
= —14.87 and 
V = 13.05. 
b. = +14.87/22.28 = +0.667. 
V(b.) = 0.1886/22.28 = 0.00846. 


Standard error = 0.092. 
As a matter of interest, the variance of b, was also determined. This 
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variance is given by 


0.1886 
[>> 27 

where the w™ are the elements of the inverse of the weight matrix; in 
other words, they are the elements of the covariance matrix of the 
y;’s. Here @’ is the unweighted mean of the z/ . 

The variance of b, was found to be 0.0310, giving a standard error 
of 0.176. 

Thus the efficiency of this estimate is 


0.00846 
0.0310 27 per cent. 


With such a large standard error, any estimate of a quantity lying 
between 0 and 1 is of little value. 
To test departure from regression, we have 


x72) = 5.302(V — U?/T) = 5.302 X 3.125 = 16.57. 


Since this value exceeds the 1 per cent point of the x’ distribution, 
the data depart significantly from the assumed linear relationship. 


6. Allowance for departure from empirical relationship 


Departure from linearity (i.e., from the empirical law of Fairfield 
Smith) may cause concern in some examples. In such cases, provided 
cost data are available, optimum plot size may be estimated with 
reasonable accuracy without the assumption that the empirical law 
holds. 


Suppose the cost of r replications is 
r(K, + 


where K, is the cost of a plot (regardless of size), K, is the cost per unit 
of plot and x is the number of units per plot. 

We then require to minimize V,/r, subject to the condition that 
r(K, + K,2) be fixed. This is equivalent to minimizing 


F(z) = (K, + K,2)V, 


with respect to x. If (x) can be determined from experimental data 


for a few values of z, its minimum may be fairly easily determined 
graphically. 


Example 2 


Johnson and Hixon [1952] have reported a 100 per cent cruise of 40 
acres of old-growth Douglas fir timber in Oregon. The data consist of 
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timber volume on each of 1600 1/40-acre plots in a 40 X 40 square. 
The analysis of variance, with stratification to eliminate systematic 
variation between sets of 8 rows and sets of 8 columns, has been worked 
out in the manner shown in the table below: 


Degrees of 
Source Freedom Mean Square 


Among 1.6-acre plots 16 277106 = Vi 
Among 0.4-acre plots in 1.6-acre plots 75 93947 = V2 
Among 0.1-acre plots in 0.4-acre plots 300 73012 = V; 
Among 0.025-acre plots in 0.1-acre plots 1200 100744 = V, 


The large mean square among 0.025-acre plots indicates competitive 
effects. Since the average diameter of trees measured was 45 inches, 
such a result is hardly surprising. 

Here 


Number of 0.025-acre 
units per plot 


277106 


When the values of the V’ are adjusted to a unit basis and plotted, 
departures from linearity appear serious. 

Johnson and Hixon also estimate the number of plots of different 
sizes which can be measured in a four-hour cruise of 40 acres. Con- 
verting these data to mean minutes per plot for plots of different sizes, 
we obtain 


| 
| Number of Mean Minutes 
Kind of Plot | 0.025-acre Units (x) per Plot (t) 
} X 1 chain | 2 7.50 
} X 2chains | 4 11.10 
} X4chains | 8 16.14 
+ X 6 chains | 12 22.16 


Vi = 64 
Vi = 138349 16 
: V3 = 89223 4 
Vi = 97869 | | 


PLOT SIZE AND CROP YIELDS 221 


Assuming cost per plot (in minutes) to be linearly related to size 
of plot, we may write 


T = K, + K:z, 


where 7’ is total cost per plot (measured in minutes per plot), 
K, is a constant (measured in minutes per plot), 
K, is a constant (measured in minutes per unit area), 
x is the number of 0.025-acre units per plot. 
From the data given above we obtain 


T = 4.9 + 1.432. 
Thus we may compute F(z): 


x (number of 
0.025-acre units) Ki + Kez V'/ze = Vz F(z) = Vz (Ki + Kez) 


1 6.33 97869 619500 
4 10.62 22306 236900 
16 27.78 8647 240200 


64 


Plotting F(x) as a function of x we find that its minimum occurs between 
x = 4and xz = 16 units. It is suggested that in this region departures 
from linearity in the relation 


log V, = log — logz 
will not be serious. Hence £ is given approximately by 
— (log Vis — log V4)/(log 16 — log 4) 
= — (log 8647 — log 22306)/(0.602) = + 0.68. 
Optimum size of plot is then 
Zo = bK,/(1 — b)K2 = (.68)(4.9)/(.32)(1.43) = 7.3 units. 


More accurate estimates of b in this region might be obtained by calcu- 
lating mean squares for plots of size 5, 8, 10, and 12 units and computing 
b using the weighting procedure suggested above. 


REFERENCES 


Johnson, F. A., and Hixon, H. J. [1952]. The most efficient size and shape of plot 
to use for cruising in old-growth Douglas-Fir timber. Jour. of Forestry 50: 17-20. 


‘te A 
| 
| 96.42 4330 417500 
. 


222 BIOMETRICS, JUNE 1958 | 


Koch, E. J., and Rigney, J. A. [1951]. A method of estimating optimum plot size 
from experimental data. Agronomy Journal 43: 17-21. 


Rao, C. R. [1952]. Advanced Statistical Metiod in Biometric Research. John Wiley 
& Sons, New York. 


Smith, H. F. [1938]. An empirical law describing heterogeneity in the yields of 
agricultural crops. J. Agric. Sci. 28: 1-23. 


ESTIMATION OF RELATIVE FREQUENCIES OF FOUR 
SPERM TYPES IN DROSOPHILA MELANOGASTER 


Marvin A. KasTENBAUM 


Mathematics Panel, Oak Ridge National Laboratory 
Oak Ridge, Tennessee, U.S.A. 


1. Introduction 


Male Drosophila Melanogaster carrying the translocation T(1, 4)B* 
[Translocation (1, 4) Bar of Stone] between the X chromosome and the 
minute fourth chromosome, also carry a Y chromosome and a normal 
fourth chromosome. Such males have been observed to produce sperm 
which carry either the proximal end of the translocation or the Y, but 
not both, and either the distal end of the translocation or the fourth 
chromosome, but not both. Thus, if A is the proximal part of the 
translocation and A’, its homolog, is the Y chromosome, and if B is 
the distal part of the translocation and B’, its homolog, is the fourth 
chromosome, then the four sperm types which may be produced are 
AB, A'B’, AB’, 

When such males are mated to attached-X females with a Y chromo- 
some, AB, A’B’, and AB’ are recoverable, and A’B is lethal. If the Y 
chromosome of the female is replaced by the proximal segment of the 
translocation, the recoverable sperm types are A’B’, A’B, and AB’, 
with AB semilethal or lethal. Results of experiments of this type 
have been reported by Novitski and Sandler [1]. 

The present paper is concerned with the derivation of estimates of 
the relative frequencies of the four sperm types as well as the variances 
and covariances of these estimates. In Section 5, the hypothesis of 
homogeneity in a specific problem is tested by means of chi-square, 
and the results are compared with those of Watson [2]. 


2. Estimation of Relative Frequencies by Method of Maximum Likelihood 


If the true probabilities of observing the four sperm types are p, 
for AB, p. for A’B’, ps for A’B and p, for AB’, where 


4 


O<p<)) (1) 
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then, in matings of the first type, where A’B cannot be observed, the 
probabilities associated with the three recoverable sperm types may be 
specified as 


for AB; q@ for A’B’; and 
(2) 
= for AB’. 
1 Ps 


Similarly, in matings of the second type, the probabilities may be 
specified as 


for A’B’; r= 3 for A’B; and 
1 Pi 


= 
(3) 


Then the probability of observing a, sperm of type AB, a, of type 
A’B’, and a, of type AB’ out of a sample of size N, resulting from the 
first kind of mating, and, simultaneously, b, sperm of type A’B’, b; of 
type A’B, and b, of type AB’ out of a sample of size N, resulting from 
the second kind of mating, where 


4+a+4=Ni, GHtetu=l, (4) 
bo +b; +h, = N2, and tm = 1, (5) 
may be written as the product of two multinomial distributions; namely 


a, Pi P2 Ps (1 Pi P2 Ds) 


6) 


It can be shown that the maximum likelihood estimates of the 
parameters of this distribution are: 


a,(b, + bs) 

+ as) + Ni(b. + 
bs(a2 + as) 

Ps + ay) + + 8) 


— (a, + b)(1 — p, — ps) 
Ps (dy + by + ay + by) 


(a, + by) | (Ni — a,)(b, + ] d 
+ a.) + Nib, ™ 


(9) 


~ (aa + be + a, + by) 


| 
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ps = 1— fi — — Ds 
= (a, + f (Ni — + Da) |. (10) 
(a, + by + ay + bx) Ld3(@2 + as) + + bu) 


3. Derivation of Asymptotic Variances and Covariances 


The matrix of asymptotic variances and covariances of the estimates 
of the parameters of the joint multinomial distribution (6) may be 
found by inverting the matrix 


A= |- (11) 


where £ is the operator of mathematical expectation, and L the logarithm 
and p; the parameters of (6), for 7 and 7’ = 1, 2, 3. 


Since E(a,) = — for j = 1,2,4, and (12) 
E(b,) = for k = 2,3, 4, (13) 
1 


the matrix A may be written as follows: 
(l1+d,) 1 1 
A=t 1 d, 1 ; (14) 


1 1 (1+ 
where 


[N,(1 — pi) + N2(1 — ps)] 


d, = (1 — Pi — Do — ps)[Ni(1 = p:)° Ps)] 


pill — p,)[N, (1 — p,) + N.(1 — ps)] 
dy = ps) (N Ds) — Nip3(1 — 


ps(1 — ps)[Ni(1 — pi) + N2( i — ps)] 
It can be demonstrated that A-', the inverse of A, is of the form 
Cu Cig Chg 


C31 C32 C33 
= Var 


pil — = 2 
~ NN Ps) [N2(1 Ds) + Nipips), 
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= Var 


D2 — — — pr — po — prs) 
(l— Ps) — + NA1 = Ps) 


+ PolN + N.pi(l — 


Var ps 


NNAL — — Ds) — pi)” + Nopips); 


= = Cov 


— P: — Ds) Ps) poli 


= Cov (fp, , ps) 


= — — p,) — + — and 


= Cov (ps Pa) 


Pops(1 Ds) 2 
N, 1- 1 N 1 — 
— p, — py — — — 


Estimates of these asymptotic variances and covariances may be found 
by replacing p, , pz. , and p; by their respective maximum likelihood 


_estimates. 


4. Application 
The data in Table 1 are results of experiments performed at the Oak 
Ridge National Laboratory and have been presented elsewhere [1]: 


TABLE 1 


NUMBERS OF PROGENY FROM MATING TRANSLOCATION-BEARING MALES TO 
ATTACHED-X FEMALES 


Male Sperm Types 
Total 


AB | A'B’ | A'B | AB’ 


1,413 | 1,029 | lethal 2,240 | 4,682 


lethal 548 346 1,287 | 2,181 
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When applied to these data, equations (7) to (10) yield 
p, = .2667, pe = .1906, ps = .1163, and ff, = .4264. (16) 
The standard errors are estimated as follows: 
S.E. = .65 X 10°’; S.E. 10°?; 
S.E. (f.) = .45 X 10°’; S.E. (p,) = 61 X 10°’. 
5. Test of Goodness of Fit 


On the null hypothesis, the observed cell frequencies represent two 
samples of size N, and N, , respectively, from the same multinomial 
population with probabilities p; , (i = 1, 2, 3, 4), where >-*_, p; = 1. 
To test this hypothesis we employ the test statistic 


(17) 


(observed — expected)* 
expected 
which is distributed approximately as chi-square, with degrees of freedom 
equal to the number cells, minus the number of constraints, minus the 
number of independent parameters estimated from the data. 
From equations (12) and (13), we may write 
E(a,) = Nipi/(l — ps) = a ; 
E(bs) = Nops/(1 — pi) = bs 
(Ni — + bz) 
= 1- = : 
E(a2) Nip2/( ps) (a, + bs + as + bs) 
A (N2 — b3)(az + bz) (18) 
= 
Nap2/( (dz + by + ag + by)’ 


1 — bs 
E(as) = — ps) = be 


2 — b;)(as 4 
E(bs) = Nops/Q — pi) = 


Thus for the data in Table 1, we find that 


E(b:) 


x’ = (18.966) + 3958.966 566.966 (19) 
1 
+ | 


with (6 — 2 — 3) = 1 degree of freedom. The probability of obtaining 
a value of chi-square as large as this or larger, under the null hypothesis, 
is approximately 23.2 per cent. We have no reason, therefore, to reject 
the null hypothesis. 
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It is interesting to note that the arithmetic results achieved in this 
section are identical with those which we would have gotten if we had 
used the techniques outlined by Watson [2]. In essence we are dealing 
with a 2 X 4 contingency table with two missing values. The two 
missing values may be estimated by a simple iterative procedure which 
Watson describes. After only a few iterations, estimates of these values 
converge to approximately 


a; = 616, and b, = 793. 


When these estimates are added to Table 1, the data are of the form 
described in Table 2: 


TABLE 2 


AB A'B’ A’B AB’ Total 


I 1,413 | 1,029 | (616) | 2,240 | (5,298) 


II (793) 548 | 346 1,287 | (2,974) 


Total (2,206) | 1,577 | (962) | 3,527 | (8,272) 


The resulting chi-square test with degrees of freedom equal to (r — 1) 
(c — 1) — 2 = 1, is identical with (19). Also the estimates of the 
parameters p,; , based on the data in Table 2, are approximately equal 
to the results in (16). It is also worth noting that the maximum likeli- 
hood estimates (16) approximately satisfy the three simultaneous 
equations in p, , 2 , and ps which result from minimizing chi-square. 
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A GENERALIZED CLASS OF CONTAGIOUS DISTRIBUTIONS* 


JoHN GURLAND 
Towa State College, Ames, Iowa, U.S.A. 


1. Introduction 


Frequently the data arising in studies of entomology and bacteriology 
cannot be described by the usual distribution functions but rather by 
some type of contagious distribution. As Feller [2] has pointed out 
there are essentially two kinds of contagion involved in such distributions. 
One type, “‘true contagion,” is due to the fact that each “favorable” 
event increases (or decreases) the probability of succeeding favorable 
events. The other type, ‘apparent contagion,” is due to an inhomo- 
geneity of 'the population. Some distributions, such as the negative 
binomial, can apparently be interpreted on the basis of both types of 
contagion. 

In a paper by Neyman [3] a class of contagious distributions is 
derived from a certain biological model which takes into account the 
fact that the distribution of larvae over the plots of a field depends upon 
the fact that the larvae are hatched from egg-masses which appear at 
random over the field. This class of distributions has been successful 
in accounting for the distribution of some insect populations (cf. Beall 
[21]). The negative binomial has also been found useful as a possible 
underlying distribution for insect populations (cf. Bliss and Fisher [23]). 
Contagious distributions have also been used in the study of accident 
and medical statistics (ef. Dubourdieu [24], Greenwood and Yule [14], 
Lundberg [25], Newbold [26].) 

A family of contagious distributions which provides a generalization 
of Neyman’s types A, B, C, was recently presented by Beall and Rescia 
[1]. The p.g.f. (probability generating function) was shown to be 


Ge) = {=m} ep {mre +) 


which for 6 = 0, 1, 2, reduces respectively to the p.g.f. of Neyman’s 
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Types A, B, C. This family of distributions was fitted by B-R (abbre- 
viation for Beall and Rescia [1]) to some biological data and found to 
afford an improvement in some cases. 

A wider generalization as shown below is afforded by the family with 
p.g.f. 


H(z) = exp {—m,} exp [mF ifa, a + B, m(z — 1)}] (2) 


where ,F,(w, v, xz) is the confluent hypergeometric function, a well- 
known series representation of which is 


v, 2) 


viv + 1) 2! vv + D@+2) 3!° 
_ylwin TO 2. 
Tw) 


It can be seen that H(z) becomes identical with G(z) when a = 1. 
Furthermore, recurrence relations analogous to those of B-R can be 
established by appealing to properties of the confluent hypergeometric 
function. 

The purpose of this paper is twofold: first, to present some generalized 
families, including (2), by means of simple considerations involving the 
appropriate p.g.f.’s; second, to consider some limiting cases of (2). In’ 
particular, it will be shown that the limiting case 8 = © in (1), con- 
sidered, inter alia, by B-R actually corresponds to a Polya-Aeppli 
distribution. As a consequence, a simpler formula than that of B-R 
is available for computing the probabilities. 

A discussion of hypergeometric functions and other mathematical 
considerations ofa technical nature are given in the appendices at the 
end of this paper. 


‘ 2 3 
w(w + 1)(w + 2) 


2. A generalized family 


Feller [2] has shown that Neyman’s family of contagious distribu- 
tions can be obtained as a generalized Poisson distribution, where the 
distribution through which the Poisson is generalized is itself obtained 
by compounding certain distributions. More specifically, in terms of 
larvae and egg-masses in a field, suppose the field is divided into plots 
of equal areas. On the assumption the larvae may come from various 
egg-masses, the probability exactly v egg-masses are represented on a 


plot under observation is supposedly given by a Poisson distribution 
with parameter m, , 


(3) 
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Now suppose the probability there are exactly n survivors is the 
same for all egg-masses and is given by a Poisson distribution with 
parameter say: 


e*—. 
n! 


(4) 


To recall briefly the arguments in [2] and [3] suppose in a particular 
egg-mass there are n survivors; then the probability s of them will be 
found on the plot under observation will be assumed to be given by the 
binomial distribution with parameters n, p say: 


(")ra (5) 


Application of a theorem by Feller [12, page 233] yields the p.g.f. 
of the distribution of larvae as that of a generalized Poisson given by 


exp [m,{g(@) — 1}] (6) 
where g(z) is the p.g.f. of the binomial distribution (5) compounded 
with a Poisson distribution (4) on the parameter n. In this case the 


parameter p in (5) is regarded as fixed for all egg-masses; hence g(z) 
is given by 


exp 


and (6) becomes 
exp [m,{exp m(z — 1) — 1}] (8) 


when m, = Ap. This is, of course, the p.g.f. of a Neyman Type A 
distribution. 

A generalization of the above may be obtained by regarding the 
parameter p in (5) as a random variable having a Beta distribution 
with parameters a, 8 (cf. Cramer [13], page 244]) 


1 B-1 
Be, a> 0, B> 0. (9) 


The result of thus compounding the distribution represented by (7) 
through values of p yields the p.g.f. g,(z), say: 


gi(z) = 3) [ — de (10) 


= {a, a + B, — 1)}. 
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(Some elementary aspects of hypergeometric functions are indicated 
in the appendices. ) 

Application of the aforementioned theorem in Feller [2] now yields 
as the p.g.f. of the distribution of larvae 


H(z) = exp [m{gi(z) — 1}] (11) 
which is the generalization given in (2), with \ replacing m, . This is 
a four-parameter family of distributions with parameters 

m, = mean density of egg-masses 
\ = mean survival number in egg-masses. 


The parameters a, 8 refer to the Beta distribution on the parameter p 
in the binomial distribution (5) giving the conditional probability that s 
of the n survivors of an egg-mass will be found on the plot under observa- 
tion. 


3. Recurrence relations 


Let P, denote the probability of an observation of r events. In (2) 
set 


iF fa,a + B, — 1)} = (12) 
Then 
H(z) = exp {—m,! exp {m,f(2)} 
= Pz’ (13) 


and the probabilities P, are obtained from the derivatives of H(z) at 
z=0. Thus 


(r) 
p 
r! 
It is obvious, of course, that 
Po = exp {—m,} exp {mf(0)}. (14) 


It can be shown (refer to Appendix A.1) that 
m, 
= r+1 (15) 
where 


F, = — (16) 
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(The notation F, should not be confused with ,/’, or F used for the 
confluent hypergeometric and the hypergeometric functions respec- 
tively.) 

From the formula for the derivative of a confluent hypergeometric 
function (formula (40) in A.2) it follows that 


= a(a + 1) ---(@+h) 
k! (17) 


k + lat+B+ k+ — M2). 


Thus, in computing probabilities P, from (15) it is possible to make 
use of tables of the confluent hypergeometric function (e.g. [4], [5], 
[6], [7], [8]) or related tables of the Incomplete Gamma Function [19| 
for obtaining the values of F, . To make use of such tables one must, 
of course, remove the negative argument —m, in (17). This is easily 
accomplished through Kummer’s formula [(41) in A.2]. 

In general, some interpolation in the tables will be necessary to 
obtain the values of F, to the degree of accuracy necessary. An example 
illustrating how the confluent hypergeometric function may be used in 
computing the probabilities is presented in Section 4. 

It is also possible to compute the values of Ff, from the following 
recurrence formula 

at+k-1 
a proof of which formula is given in A.4. On setting a = 1 it is evident 
that (18) reduces to B-R’s formula (24). [There is a slight error in 
their formula (24). The coefficient of F,-. should be m./(k — 1) 
instead of m,/k. In their computations in Table XI, loc. cit., they 
use the proper coefficient m./(k — 1).] 

As formula (18) suffers from a defect similar to that mentioned in 
B-R, namely, that rounding errors in the earlier values of /, may 
seriously affect later values of F, , it is desirable to compute the values 
of F, directly and to use formula (18) as a check on computed values. 


4. Example 


For convenience we consider the insect data of B-R’s Table II. 
For the case a = 1, 8 = 2 the estimates of the parameters obtained by 
the method of moments are given as 


Mm, 1.114995069 


m, = 8.023738236. 
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They suggest this large number of decimal places to compensate 
for the round-off errors in their computations, especially those due to 
the use of formula (18) with a = 1, 8 = 2. 

To compute the values of F, we shall use the alternative method 
based on formula (17) which involves explicitly the confluent hyper- 
geometric function and which is much less sensitive to round-off errors. 
Specifically, for the case a = 1, 8 = 2, and the estimated value of m, 
given above, it is required to determine the values of 


.F\(2 + k, 4 + k, —m,) 


for successive integral values of k from 0 through 14. 
Application of Kummer’s formula [(41) in A.2] yields 


iF (2 + k, 4+ k, — Mz) = F,(2, 4+ k, M2) (19) 


which is advantageous in that the first of the three variables in the 
confluent hypergeometric function appearing on the right-hand side 
remains constant for the different values of k, and further, the third 
variable is positive. 

For the values of the variables w, v, x in ,F,(w, v, x) which are re- 
quired in this example it is more convenient to use the Tables of the 
Incomplete Gamma Function [19] in computing ,F,(w, v, x) than to 
interpolate in and extend the tables of the confluent hypergeometric 
function included in the list of references. In fact, until such time as 
more extensive tables of the confluent hypergeometric function become 
available, it is probably preferable in most cases to use the Tables of 
the Incomplete Gamma Function in such computations, especially 
for low integral values of the variable w. 

Adaptation of formula (43) in Appendix A.2 enables us to compute 
the values of ,F,(1, k + 2, m.) from the following: 


iF,(1,4 + k, = (8 + at (20) 


where the estimated value of m, to be used has been indicated above. For 
such a value of m, interpolation is necessary in the Tables of the Incom- 
plete Gamma Function. Bessel’s interpolation formula (ef. [22]) 
involving differences up to the fourth order has been found adequate, 
and the Tables [19] are very convenient for this purpose since the second 
and fourth order differences are tabulated. 

From the values of ,F,(1, 4 + k, m) it is a simple matter to compute 
the values of ,F,(2, 4 + k, m.) through application of the recurrent 
formula (42) in Appendix A.2, and noting that ,F,(0, v, x) = 1. 

In Table 1 column (1) are shown the values of F, for the case a = 1, 


af 
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TABLE 1 
(a = 1,8 = 2) 
(1) (2)* (3)** 

k m Fi Fy 
0 . 1872316 . 208762368 . 1872316 
1 .3131714 . 349184607 .3131714 
2 . 3804478 .424197564 3804479 
3 . 3952149 . 440663069 3952152 
4 . 3680596 410385813 . 3680606 
5 .3131510 . 349164527 3131534 
6 . 2456656 . 273921106 2456702 
7 . 1786759 . 199231896 1786841 
8 . 1209571 . 134881271 1209703 
9 .0764668 .085282038 0764865 
10 .0452780 .050516297 0453063 
11 .0251840 .028121143 0252209 
12 .0131928 .014762919 0132403 
13 .0065261 .007342910 .0065856 
14 .0030559 .003488466 .0031287 


*The values in column (2) are reproduced from B-R’s Table XI. 


**The values in column (3) are obtained from those in column (2) by dividing by m: = 1.114995069. 


TABLE 2 
= 1,8 = 2) 
(1)* (2)** 

P, 

0 .41824 41824 

1 .08731 08731 

2 .08214 08214 

3 07502 07502 

4 06642 06642 

5 05700 05700 

6 .04750 .04750 

.03854 .03855 

8 03060 .03060 

9 02387 02388 
10 01839 01839 
ll 01384 01403 
12 01062 01063 
13 00799 00802 
14 00597 00601 


*The values in this column are computed from the values of Fy in column (1) of Table 1. 
**The values in this column are reproduced from B-R’s Table XIII. 
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8 = 2, computed through the values of the confluent hypergeometric 
function as described above. These values are correct to within one 
unit in the seventh decimal place. In column (2) of Table 1 are shown 
the values of m,F, which appear in B-R’s Table XI, and in column (3) 
are shown the values of F, obtained through dividing the numbers in 
column (2) by m, = 1.114995069. It is apparent that the magnitudes 
of the errors in the values of F, in column (3) increase substantially 
with increasing values of k. 

To compute the probabilities P, for the case a = 1, 8 = 2 we apply 
the recurrence formula (15) to the successive values of F, in column (1) 
of Table 1. These values are listed in column (1) of Table 2. For 
comparison, the values of P, given in B-R’s Table XIII are also listed, 
and it is seen that the agreement is good. This would indicate that the 
round-off errors through using formula (18) do not have a serious effect 
on the computed values of P, . In computational work one may risk 
using the recurrent formula (18), but as B-R point out, it is necessary 
to retain a considerable number of decimal places for internal accuracy. 

For the purpose of illustrating how to fit members of the family (2) 
with a ~ 1, the computations for two other cases employing the same 
data as above are summarized here. These cases are a = 2, 8B = 3 
and a = 3,6 = 4. 

In Table 3 the expected and observed frequencies are shown for the 
case a = 2, B = 3, along with the computed values of the terms F, 
and the probabilities P, . These values have been computed by the 
same method as that described above for a = 1, 8 = 2, and by making 
use of recurrence formulae for the confluent hypergeometric function. 
In Table 4 the corresponding entries, similarly computed are given for 
the case a = 3, 8 = 4. In Tables 3 and 4 most of the values of F, are 
correct to within a small error in the seventh decimal place, except for 
the last few values in Table 4. The values of F, computed by means of 
the recurrence formula (18) have been inserted in parentheses in Tables 
3 and 4. The agreement between the two methods is comparable to 
that for the case a = 1, 8 = 2 presented in Table 1. Also, as for the 
case a = 1, 8 = 2 the probabilities and expected frequencies computed 
on the basis of these approximate F, values are practically the same as 
those appearing in the table and based on the more exact F, values. 
The values of x’ based on grouping similar to that in B-R’s Table II 
are exhibited along with the corresponding probabilities of exceeding 
these values of x’. 

The problem of how to obtain the most efficient estimators of the 
parameters in a distribution of the form (2) is beyond the scope of the 
present paper. For the purpose of illustration of the techniques in- 
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TABLE 3 
(a = 2,8 = 3) 


Expected Observed 
F, Py frequency frequency 


0 . 1632712 (.1632712) .43771 24.5 24 
3479885 (.3479885) .06640 3.7 
2 .4749821 (.4749820) .07580 4.2 4 
3 5172230 (.5172229) .07538 4.2 5 
4 . 4839899 (. 4839896) .06890 3.9 3 
5 4028450 (.4028442) 05941 3.3 2 
6 3040002 (3039983) 04912 2.8 1 
7 .2104938 (.2104898) .03941 2.2 1\ 
8 . 1348416 (. 1348343) 03098 2 
0804104 (.0803975) 02401 1.3 1 
10 0448570 (0448362) 01843 1.0 2 
0235044 (.0234719) 01404 0.8 3 
12 0116080 (.0115599) .01061 0.6 1 
13+ .02980 1.7 1 


x? = 5.33, P(x?) = .50. 


TABLE 4 
(a = 3,8 = 4) 


Expected Observed 


k P, frequency frequency 
0 . 14380719 (.1430719) .43943 24.6 24 
1 . 3509552 (.3509552) .05452 3.0 6\ 
2 .5158642 (.5158642) .07025 3.9 4) 
3 . 5806688 (.5806683) .07397 4.1 5 
4 . 5468752 (.5468739) .06906 3.9 3 
5 .4500262 (.4500214) .05967 3.3 2 
6 .3317386 (.3317257) .04904 2.7 1 
7 . 2225779 (.2225477) .03906 2.2 1 
8 . 1374319 (. 1373666) .03053 2 
9 .0787297 (.0784477) .02149 1.2 1 
10 .0421119 (.0415015) .01807 1.0 2 
11 .0211226 (.0200224) .01372 0.8 3 
12 .0108483 (.0081036) .01033 0.6 1 
j 2.8 1 


x? = 6.84, P(x?) = .34. 


volved, the values of a and 8 have been assigned (cf. B-R) and the 
values of m, , m, estimated by the method of moments. The problem 
of estimation in multi-parameter families will be presented in a later 
paper. It might be remarked here, however, that on the basis of certain 
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empirical studies it has been found that the first two sample moments 
in conjunction with the sample frequencies of low counts show possi- 
bilities of effective estimation of the parameters. 


5. Some limiting cases 


The form of the generalized family (2) will now be considered for 
limiting values of the parameters a, 8. 
(a) B— &, a fixed 

Now 


and 


_ 1 
@+Ae+64+ 1) 
Hence, keeping a fixed 


Ta@+r) T@+8) 

= 

T@ 
and ,F,{a, a + B, m.(z — 1)} — 1 which is the p.g.f. of a degenerate 
distribution corresponding to the constant zero. 

It is interesting to note that as the family (2) approaches this limit 
it passes through a family of generalized Polya-Aeppli distributions 
for large values of 8. This can be seen as follows: 

Since 
T@+8+7) 

+ 8) 


it follows, for a fixed, 


T(a + 8) ‘|- 


Hence, for large 8, we can write approximately 


= exp | m{(1 mie > (21) 


lim 


Bow 


7 
| 
I 
} 
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iF; {a,a + B, 2 I'(a) 
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which is the p.g.f. of a generalized Polya-Aeppli distribution (cf. Skellam 
[9]). When @ = 1 this is the p.g.f. of an ordinary Polya-Aeppli distri- 
bution. 


(b) a— B fixed 


From the expressions for '(a + r)/T'(a) and + 
in (a) above it is easy to see that for 6 fixed 
Tet+f 


Consequently 
iF {a,a + B, m(z — 1)} exp {m,(z — 1)} 
and H(z) approaches the p.g.f. of the Neyman Type A distribution 
exp [m,{exp m.(z — 1) — 1}]. (22) 
(c) a> ©, ~, B/a = 
According to this case the values of a and 8 are allowed to approach 


infinity under the restriction B/a is a constant equal to y. 
From the expressions in (a) above it can be seen that 


T@) 


B/a=¥ 


Hence ,F,{a, a + 8, m.(z — 1)} — exp {m, (z — 1)/1 + y} 
and 


H(z) exp exp (zg — 7 (23) 
the p.g.f. of a Neyman Type A distribution. 
(d) with Mi » He fixed 


One of the apparent reasons for this restriction (cf. B-R [1], Ans- 
combe [10]) is the use of the first two sample moments in the estimation 
of the unknown parameters. 

From (2) it is seen that the factorial cumulant generating function 
(cf. Aitken [17, page 23]) is 


aa +1) (ma)* 
a+ 2! 
Hence, the first two factorial cumulants are 


log + u) = m 


mu + 


My + 1) 


= 


& 
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q 
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and the ordinary cumulants are 


Define 


c= 


and it is now possible to write 


a atBt+le’ 


By an argument similar to that in (a) above it follows that 


Ta@a+r) T@+t+8) 
Ta) 


— = fe — 


fixed 


Consequently 


Fifa, a+ 8, 0} > exp tole 


and 


H(@) — exp“ [exp {e@ — 1)} 11, 
the p.g.f. of a Neyman Type A distribution. 
(e) B— @, with a, py , we , fixed 


A similar argument to that used in the preceding cases yields 


Consequently 


c 
iF {a,a@ + B, mz — 1)} - 


and 


| 

which is the p.g.f. of a generalized Polya-Aeppli distribution. Note the 
similarity of this result with that in (a) above, equation (21). 
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The case a = | corresponds to the limiting case considered in B-R. 
For this case the limiting distribution is Polya-Aeppli, and this result 
will be utilized in Section 6. 


(f) a— ©, with B, yf , ue fixed 


Similar arguments to the above yield 


T@+r) T@+8) 
ao T@) 


Consequently 


{m(z — = — 1)}’. 


iF, {a, a + B, m(z — 1)} > exp {ee — 1)} 


and 
H(z) — exp {exp — 1) — 1}, (26) 


which is the p.g.f. of a Neyman Type A distribution. Since the value 
of 8 does not enter into the limiting distribution it is not at all surprising 
this limiting distribution is the same as the one obtained in (d) above. 

A summary of all the limiting forms considered in (a) through (f) 
is given in Table 5. 


) TABLE 5 
LimitinG Form or DistRIBUTIONS IN GENERALIZED (2) 
Restriction Approach Limiting form of distribution 
a= oo Neyman Type A 
B = 
B = © Degenerate 
B large Generalized Polya-Aeppli 
a= o Neyman Type A 
a= o Neyman Type A 
B= 
a B = @ Generalized Polya-Aeppli 
Neyman Type A 
5) 6. An example 
As evident from part (e) of Section 5 the limiting case a = 1,8 = @ 
he considered by B-R under the restriction the first two moments are 


fixed corresponds to a Polya-Aeppli distribution with p.g.f. 


- 
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exp 7H! 5@ i]. (27) 


As in Section 3, write 
-1 
{@) = {I »} 


and 


as is evident from (14) and (15). 

There is, however, a simpler recurrence formula than (28) for the 
probabilities in a Polya-Aeppli distribution (ef. Evans [11]). This can 
be seen as follows: 

From the above _— for F, it follows that 


= 


— Pat Sh, =0 


where 


Combining this with (28) yields 


c(r + 2) e(r + 1) a 1) 
4Au; 


+= 


2A 


[FP — — FP, . 


Upon simplification this becomes 


_ 3A Qui 


Formula (29) was applied to the data of several tables in B-R and 


(29) 
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the agreement with their computed values of expected frequencies for 
the limiting case was extremely close. As an example, the data of their 
Table VII is given below in Table 6, together with the expected fre- 
quencies using both methods, and the agreement is seen to be very 


TABLE 6 
ComPuUTED FREQUENCIES FoR B-R Taste VII 


Observed Expected frequency 
frequency 


Expected frequency 
using (29) 


6 32 33.0 33.0 
7 17 22.9 22.9 
8 4 15.6 15.6 
9 7 10.4 10.4 
10 7 6.9 6.9 
a1 2 4.6 4.5 
12 3 3.0 3.0 
13 3 1.9 1.9 
14 1 1.2 1.2 
1 8 8 
1 5 5 
2 3 3 
1 
0 


close. In using (29) it is suggested as a note of caution that the coeffi- 
cients of P,,, and P, should be small, and in any case should not exceed 
unity in absolute value. 


7. Other generalized families 


Generalized families analogous to (2) are easily obtained by similar 
considerations. Suppose, for instance, the parameter \ in (4) is no 
longer assumed constant over all egg-masses but is regarded as a random 
variable with a Gamma distribution, with probability density, say 


5° 0-1 


= To ° A> 0, 0, 6> 0. (30) 


As is well known, the distribution obtained by compounding (4) in 
this manner (cf. Greenwood and Yule [14]) is that of a negative binomial 


| a 
0 206 226.2 226.2 
2 128 113.9 113.9 ne 
3 107 87.7 87.7 ak. 
4 71 65.0 65.0 el 4 
5 36 46.9 46.9 4 4 
nd 
| 
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random variable. A simple form for the p.g.f. of such a variable, which 
exploits the analogy with the ordinary binomial variable is (cf. Wilks 
[15, page 54]) 
— (31) 
where 
m>0; 


If the binomial distribution (5) is compounded with a negative 
binomial distribution (31) on the values of n, the counterpart of (7) 
becomes the p.g.f. 


{1 — pp — (32) 


If, further, the parameter p of the binomial in (5) is regarded as a 
random variable with a Beta distribution given by (9), the counterpart 
of g:(z) in (10) becomes g.(z), say, where 


= F{k, ,a,a + 8, pz — 1)} 
where F(u, v, w, x) is the ordinary hypergeometric function, a well 
known series representation of which is 
u(u + + 1) 2” 
ww + 1) 2! 
Tw) 2. 
Teo Tv 
The counterpart of (11) is now the generalized family with p.g.f. 


exp [m, {g2(z) — 1}]. (34) 


Other types of generalized families are also possible. For instance, 
if the Poisson distribution of (3) is replaced by a negative binomial 
with p.g.f. say, (g2 — pz)~“*, and the parameter p in (5) is regarded 
as a random variable with a Beta distribution (9), the p.g.f. of the 
resulting generalized family would be 


{1 + Pe P2 iF, {a, a B, AZ (35) 


If, in the preceding case, the Poisson distribution of (4) is also 
replaced by a negative binomial with p.g.f., say, (q. — p.2)~—"*, the 
p.g.f. of the resulting generalized family would be 


[l P2— pF {ky ,a,a+ B, plz (36) 


+ 


F(u,v,w,2z) = 1 + s+ 
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For brevity and simplicity we have confined our attention to the 
family (2). Analogous results are, of course, obtainable in a similar 
manner for the other families above. 

It may also be remarked, in passing, there are other types of com- 
pound and generalized distributions based on the negative binomial, 
which have been considered elsewhere (cf. Gurland [16]). 


8. Summary 


By appealing to the notions of compound and generalized distri- 
butions various classes of contagious distributions are obtained. One 
of these classes includes as a subclass the family of distributions con- 
sidered by Beall and Rescia. Properties of the confluent hypergeometric 
function are utilized in establishing recurrence relations for the prob- 
abilities. Some limiting cases are discussed; in particular it is pointed 
out the limiting case 8 = © of Beall and Rescia is a Polya-Aeppli 
distribution with a consequent simplification of the computation of 
probabilities. The question of estimation of parameters in multi- 
parametric families of contagious distributions will be considered in a 
later paper. 
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APPENDIX 
A.1 Proof of recurrence formula (15) 
Consider the p.g.f. 
H(z) = exp {—m,} exp {m,f(2)}. 


The first derivative is given by 


H"(z) = [exp {—m,} exp {mif(2)} = 
The (r + 1)st derivative is given by 


= m, exp t-m)(4) exp m,f(2)]. 


Now apply Leibnitz’ formula (cf. Goursat-Hedrick [20, page 28]) for 


= 
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the r-th derivative of a product u(z)v(z), say 


This yields 


HO" = m, exp {=m} @lesp 


k=0 


Pra = and F, 


Consequently 


Pra 


which is the required recurrence relation. 


A.2 Some remarks on the confluent hypergeometric function 


The differential equation 
dM 
ta + =0 (37) 


is satisfied by the confluent hypergeometric function. Two linearly 
independent solutions of this equation are (cf. Copson [18]) 


iF,(a, y, x) and F,(a, 7, 2) 


where 


iFi(a, Y 2) (a) + r) r! (38) 


(The solution x'~” ,F,(a, y, x) is valid if y is not an integer.) 
The above series is absolutely and uniformly convergent for all real 


values of a, y, x except for y a negative integer or zero. 
By utilizing the definition of the Beta function, 


Bia, b) = — dt, a> 0, b> 0, 


k 
d (7) (k), (r—k) 
twee = (7 | 
But 
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it can be shown that 


iF (a, 1,2) = Ty — a)F@) Jo t) t dt (39) 


where a and y — a must be greater than zero. 


The following well known relations can be derived from all the above 
facts: 


Fila, 7,2) = + 1,4 1,2), (40) 
iF (a, 7,2) = — a, 2), Kummer’s formula, (41) 

aF (a + 1,7, 2) = + 2a — 2) 

+ (y a) = 1, Y x), 


+ 1,2) = dt. (43) 
0 


(42) 


There are several other recurrence relations and many interesting 


properties of the confluent hypergeometric function which are presented 
in Erdelyi [27]. 


A.8 Some remarks on the hypergeometric function 


The mathematical treatment of the ordinary hypergeometric 
function, F(a, 8, y, x) say, is similar in many respects to that of the 
confluent hypergeometric function considered above. 

The hypergeometric series 


2 


is a solution of the differential equation 


xe — 1) 5+ =0. (45) 


The series converges absolutely when x < 1 and diverges when 
x > 1. The function F(a, 8, y, x) can be defined for all x by analytic 
continuation in the complex plane (cf. Copson [18)]). 

From the definition of the Beta function above it can be shown that 


fory — B > 0. 
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Other interesting properties of the hypergeometric function are pre- 
sented in Erdelyi [27]. 


A.4 Proof of the recurrence relation (18) 
Since , F(a, y, x) satisfies (37), and since (40) implies 


d’,F (a, y, x) + 1) 


dx’ = vy + 1) iF (a + 2, 7 + 2, 2), (47) 
it follows that 
File + 2,7 +2, 2) pe 
(a + (48) 
_ » 
(a + Dz iF y, 2) = 0. 
Apply this relation to the function 
ala + 1) (a+k — 2) 
and we obtain 
F, k F,-, + k(k 1) Fy-2 0. (50) 
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RECENT EXPERIENCES WITH AREA SAMPLING FOR 
AGRICULTURAL STATISTICS 


R. E. Vickery 


Agricultural Marketing Service, U. S. Department 
of Agriculture, Washington, D. C., U.S.A. 


In June 1952, a Special Subcommittee of the House Agriculture 
Committee made a number of recommendations on the procedures used 
by the Crop Reporting Service of the Department of Agriculture. One 
of the most specific was that a research unit be established to conduct 
a program of experimentation on a large scale. That research unit was 
activated. A Panel of Consultants, consisting of prominent statisticians 
and agricultural economists, was also appointed to review studies 
already in progress and to help frame a long-range program. 

Several meetings with the Panel led to the conclusion that the over-all 
objectives could best be met by a continuing research program on both 
sample survey procedures and objective crop forecasting methods. 
There was general agreement that the studies should be conducted on 
a fairly extensive scale in the field, under practical operating conditions, 
to provide a better indication of how those procedures would behave 
when put to use on a large scale and also to train the nucleus of a staff 
for full operations later. 

A plan was developed whereby such a program was to be started in 
10 southern states, and expanded to the rest of the country, region by 
region, as fast as practicable, over the next few years. In June 1954, 
the first year of the expanded program, a sample of 700 area segments 
covering about 3,000 farms in 100 counties of 10 southern states was 
interviewed to obtain data on crop acreages and livestock inventories. 
Subsamples of the farms enumerated in June were visited periodically 
during the growing season for plant observations to be used in objective 
yield forecasting studies for certain crops. Mail surveys with non- 
respondent follow-ups were also made in the fall and winter to parallel 
the Division’s other acreage and production, and livestock surveys. 

In this first June Survey, the place of residence of the farm operator 
was used as the criterion of whether or not a farm was “in” a selected 
area segment. If the operator lived inside a selected segment, his farm 
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was in the sample regardless of where the land was located. Farm 
managers were considered to be operators, but croppers (tenants 
owning no work stock or tractor power) were not regarded as operators. 
Cropper operations were included as a part of the landlord’s unit to 
avoid difficulties encountered in previous surveys of this kind. 

As many operators of farms live in urban areas, it was necessary to 
draw sample segments from such areas as well as from the open country. 
It was found that area segments of the size used in the survey often 
picked up “pockets” containing extremely large numbers of farm 
operators in urban areas. 

Respondents were asked to report for all land they operated, includ- 
ing that assigned to croppers, but to exclude data for land rented out 
to bona fide cash or share renters. The survey disclosed that even though 
respondents correctly deducted land rented out in reporting the scope 
of their over-all operations, they sometimes reported crops on such 
land in reporting the individual crop acreages. Such over-reporting 
seemed to be particularly serious for crops grown on share, such as cotton 
and corn. Direct expansions of the sample data for acres of cotton 
and corn planted were 134 per cent and 122 per cent, respectively, of 
official Crop Reporting Board estimates. For cotton the official estimate 
is largely based upon measured acres reported under the acreage control 
program. Ratio estimates from the survey for these crops, derived 
from the data reported for 1953 as well as for 1954, were in line with 
official estimates. Survey estimates of all cattle were above official 
estimates, while the survey estimates of all hogs were lower. 

A discussion of the validity of official estimates is outside the domain 
of this paper. But with respect to cotton and corn acreages, there is 
good reason to believe that the official series of annual estimates are at 
about the right levels. Cotton acreages are enumerated completely 
each year under the acreage allotment program. Although such com- 
plete data on corn are not available each year, the general level of 
plantings can be judged fairly accurately from Census data. For the 
region as a whole, Census data on acreages harvested show no large up 
and down fluctuations from one census to another. For that reason, 
coupled with comparisons of 1954 survey results with those obtained 
with a better interview procedure in 1955 for the same farms, there is 
ample evidence that the 1954 survey results were too high. The 1955 
data also confirmed the hypothesis that double-reporting on these 
crops was responsible for the bias. 

Official estimates of cattle and hogs are subject to greater uncertainty 
than those on cotton and corn acreages. The differences between 
survey results and official estimates are mentioned here as a matter of 
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interest but without any implication that the official estimates neces- 
sarily provide a hard and fast standard of comparison. 


Sampling errors of the direct expansions for selected items in the 
survey were: 


TABLE 1 
SampLinG Errors FoR SELECTED ITEMs, JUNE 1954 


Sampling 
Item 


Cotton planted 
Corn planted Acres 8.9 
All cattle No. 11.3 
All hogs 


In June 1955, the survey was repeated on the same area sample as 
for 1954, except that the number of urban segments was tripled and 
their size reduced proportionately. The questionnaire was designed 
in such a way that each farm operator was permitted to report for his 
entire unit as he himself pictured it in his mind, rather than attempting 
to limit him to an arbitrary statistical definition. In the course of the 
interview, he was questioned about the various types of tenants to 
whom he rented land or who were living on his place. Information 
was obtained about the crops grown by each tenant and the type of 
tenancy established. Similar data were obtained on livestock, separately 
for the landlord and for each tenant. When the data were later tabu- 
lated in the office, it was possible to classify each individual according 
to his proper status as a farm operator and to allocate the data reported 
for the entire operation to each operator without duplications. The 
“manager” category was dropped as a separate class in the definition 
of ‘farm operator.” 

A supplementary list of 1,000 large farms was enumerated in 1955, 
in addition to the area sample, in a further attempt to reduce sampling 
errors. Problems of farm definition and identification of place of 
residence of large farm operators were encountered. These were trouble- 
some and could cause other errors that would nullify the benefits of 
reduced sampling errors. Changes in the questionnaire did eliminate 
over-reporting of crop acreages, but reduction in sampling errors from 
modifications in the sample design was less than anticipated. Sampling 
errors for selected items in the 1955 survey are compared with those in 
the 1954 survey in the following table. 
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TABLE 2 
SAMPLING ERRORS FOR SELECTED ITEMs, 1954 AND 1955 


Item Unit June 1954 June 1955 
Per cent Per cent 
Cotton planted Acres 14.3 13.7 
Corn planted Acres 8.9 7.3 
All cattle No. 11.3 8.9 
All hogs No. 13.5 9.4 


In June 1955, an entirely different approach was also tried on a small 
scale in 100 additional area segments in these same 10 southern states. 
Instead of working with the farm as the unit of observation, a “closed- 
segment” approach, which had previously been tried by the Division 
in an acreage survey in Iowa and Indiana some years ago, was invoked. 
This involves accounting for the use currently being made of all land 
falling wholly within the boundaries of every selected area segment, 
regardless of farm boundaries or the location of the operator’s residence. 
Segment boundaries were delineated on aerial photographs, data on 
crop acreages recorded for every cultivated field, and livestock numbers 
recorded, tract for tract, for livestock actually present within the segment 
boundaries. As this approach disregards the farm as the unit of observa- 
tion entirely, it eliminates the need for complex farm definitions and 
for rules specifying which farms must be considered as “‘in” the sample. 
Such definitions and rules tend to be complicated and give rise to error. 
Results from this limited study in 1955 were so encouraging that the 
entire June Survey was conducted by that method in 1956. 

Sampling errors for selected items in this 100 segment survey were as 


TABLE 3 


SAMPLING ERRORS FOR SELECTED ITEMS FROM CLOSED- 
SEGMENT SuRVEY, JUNE 1955 


Item 


Unit 


Sampling 
error 


Cotton planted 
Corn planted 
All cattle 

All hogs 


Acres 
Acres 
No. 
No. 


Per cent 
10.7 
11.9 
14.9 
17.8 
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shown in Table 3. When comparing these results with those in Table 2, 
allowance must be made for the fact that the sample was only 1/7 as 
large in this instance. 

In 1956, the area covered by the survey was extended to 11 north 
central states and to Virginia and Kentucky, in addition to retaining 
the 10 southern states covered the previous two years. The total sample 
consisted of 1,100 area segments. There was some doubt about the 
practicability of the closed-segment approach for recording livestock 
data because livestock are free to move about. For that reason, live- 
stock numbers were reported both on a closed-segment basis and on an 
open-segment (by farms) basis. In the closed-segment work, data 
are recorded only for livestock actually present within the segment 
boundaries at the time of the interview. In the farm or open-segment 
approach, livestock numbers are recorded for an entire farm if the farm 
operator lives within the segment boundaries. 

Preliminary results from this survey show no significant differences 
in estimated livestock numbers from the two approaches in the north 
central states. In the southern states, estimates from the closed-segment 
approach tend to be somewhat higher than with the farm approach. 
Part of the difference may be caused by a tendency to include too many 
animals in recording the numbers within the segment boundaries in 
the closed-segment approach. But it is also possible that there is some 
under-reporting of livestock when farmers are asked to report on all 
their holdings. 

Some idea of the effect of the closed-segment approach on sampling 
errors can be gained by comparing results from the 1955 and 1956 
surveys in the southern states. The 1955 sample was.somewhat larger 
than the 1956 sample. In addition, the 1955 sampling errors include 
benefits obtained from enumerating 1,000 large farms. 

In 1955, the segments were distributed over 100 counties, whereas 
they were spread over 325 counties in 1956. This wider dispersion, no 
doubt, reduced between-county contributions to sampling errors 
appreciably. The comparisons are shown below. 

In general, experience in these surveys indicates that the closed- 
segment approach works as well in recording livestock data as for 
land-use items. In addition to other advantages, it results in a tre- 
mendous gain in statistical efficiency for land-use items by limiting the 
data to land actually within the segment boundaries. The range within 
which data reported for individual segments can vary is greatly reduced. 
Although the reduction in sampling error varies from item to item, a 
conservative generalization is that one segment, using the closed-seg- 
ment approach, yields at least as much information as two segments 
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TABLE 4 


SaMPLiING Errors For SELECTED ITEMS IN THE SOUTHERN 
Sratss, 1955 anp 1956 


Sampling errors 
Item Unit 
1955 1956 
Per cent Per cent 

Cotton planted Acres 13.7 7.2 
Corn planted Acres 7.3 5.2 
All cattle No. 8.9 5.4 
All hogs No. 9.4 6.9 


using the farm approach. Sampling errors for livestock numbers are not 
reduced appreciably by the closed-segment approach, except for cattle. 
As livestock data were obtained on both the ‘“‘closed-segment”’ and 
“open-segment”’ (or farm headquarters) approaches in the 1956 survey, 
some indication of relative efficiency for these items can be gained by 
comparing sampling errors. The comparison is shown in Table 5. 


TABLE 5 


SAMPLING Errors oF Livestock EstTIMATES BY ‘‘CLOSED’’ AND 
“OpEN’’ SEGMENT APPROACHES 


Efficiency of 
Closed Open closed segment 
Item segment segment over open segment 
Per cent Per cent Per cent 

A—Southern States 

All cattle 5.4 6.4 140 

All hogs 6.9 6.6 91 

All sheep 14.4 9.1 40 

Hens & pullets laying age c (mek 7.6 115 
B—North Central States 

All cattle 5.4 6.9 163 

All hogs 5.8 5.7 97 

All sheep 10.1 11.7 134 

Hens & pullets laying age 5.1 5.3 108 


As might be expected there was little gain in efficiency from the 
closed-segment approach for items closely associated with farm numbers, 
such as hogs and chickens. Sheep are a rare item on farms. 
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All sampling errors quoted in this report were computed by treating 
the data as a sample of counties. The sample design involved making a 
random selection of counties within type-of-farming areas within 
states and selecting area segments at random within those counties. 
Counties were selected with probabilities proportional to numbers of 
farms. The sampling rates with which area segments were selected 
within counties were inversely proportional to the probabilities with 
which the counties themselves were selected. With such a design, 
the sampling error can be computed very simply by first computing a 
per-segment average for each county and then computing the between- 
county variability of those averages within type-of-farming strata 
within states. The sampling variance of the per-segment average for 
the entire region was obtained by dividing the between-county variance, 
computed as indicated above, by the number of counties in the survey. 

Field costs averaged about $53 per segment in 1954 and 1955 when 
the “‘open-segment” or farm headquarters approach was used. Such 
costs were about $65 per segment in 1956. 


™ 


ESTIMATION OF MISSING VALUES FOR THE ANALYSIS 
OF INCOMPLETE DATA 


G. N. WILKINSON 


Commonwealth Scientific and Industrial Research Organization 
Division of Mathematical Statistics 
Adelaide, Australia 


1. INTRODUCTION 


When only a few observations are missing from data that otherwise 
conform to a planned experimental design, the fitting of a linear model 
to the data by the principle of least squares is most simply carried out 
using estimated values for the missing observations, as in this way 
most effective use is made of the remaining symmetry in the experi- 
mental results. 

A comprehensive account of this procedure was given by Yates [12] 
in 1933. Firstly, the data are completed by inserting estimates of the 
missing values, which are determined so that the residual sum of squares 
for the completed data shall be a minimum. The correct estimates of 
treatment effects and other parameters in the linear model are then 
given by the standard formulae for the experimental design, and the 
correct residual sum of squares for the incomplete data is given by the 
standard analysis of variance (with the degrees of freedom appropriately 
reduced). However, the other components of this analysis are only 
approximations (though usually reasonably good) to the corresponding 
components of a correct least squares analysis of variance, since the 
latter involves in principle the fitting of one or more auxiliary models, 
for which appropriate estimates of the missing values would be necessary. 
It should also be noted that the usual formulae for standard errors of 
treatment comparisons need to be modified to take into account the loss 
in precision due to the missing values. 

The non-standard aspects of this missing value procedure will be 
discussed by the present writer in two papers: the present paper, dealing 
with the estimation of missing values (which enables completion of 
the first phase of analysis), and a second paper [11] which deals with 
correction of the analysis of variance and the derivation of standard 
errors. In connection with incomplete block designs, lattice squares, 
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and split-plot designs, it should be mentioned that the procedure 
discussed in these papers is strictly appropriate only for the intra- 
block and sub-plot analyses, respectively. For the recovery of interblock 
information a weighted analysis is necessary, the missing value pro- 
cedure for which has been given by Cornish [3], and will be further 
discussed by the writer in a future publication. (For a discussion of 

missing values in split-plot designs, see Anderson [1].) 

The main emphasis of the present paper is on the following points: 

(i) Equations for missing values can be simply determined. The 
currently recommended procedure, due to Yates [12], for determining 
the estimates of several missing values, is by the iterative application 
of the appropriate formula for a single missing value. Yates noted 
that, in principle, equations for the missing values could be derived 
by differentiating the standard residual sum of squares with respect 
to unknowns representing the missing values, but apparently did not 
observe that the equations could be determined in a more direct way: 
Each unknown is simply equated to its estimated value (a function of 
the unknowns) derived from the formally completed data. The latter 
procedure was mentioned by DeLury [4], but otherwise has not received 
much attention in the literature. 

(ii) Solution of the equations by matrix inversion is usually relatively 
simple. The main reason for considering solution by this method, is 
that the inverse of the matrix of coefficients is also required in order 
to determine correct standard errors for treatment comparisons. The 
necessary formula was given by Tocher [9] and will be further con- 
sidered in [11]. Apart from this important aspect, however, it will be 
found that in the usual situation, in which the missing values are few 
in number and have a reasonably simple configuration within the design, 
equations for the missing values can be determined and solved with less 
computation than would be required for the iterative procedure referred 
to earlier. 

The considerations below fall naturally into three parts and are 
dealt with in the following sections: 

2. THEORETICAL, giving mainly a derivation of the missing value 
equations. 

3. SETTING UP THE MISSING VALUE EQUATIONS, in which 
methods of computation are discussed. 

4. SOLUTION BY MATRIX INVERSION, in which some well 
known matrix inversion formulae are given and their application 
discussed. The solution of singular equations is also discussed. 
Section 2 also deals with the missing value procedure for an analysis 

of covariance. Detailed illustration of the procedure in this situation 

will be found in an earlier paper by the writer [10]. 
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2. THEORETICAL 
Assumptions 


The present paper is concerned with the formal application, rather 
than the statistical appropriateness, of the method of least squares 
(unweighted). Accordingly we shall assume a linear model for the 
data, in which all relevant experimental factors are formally represented 
by parameters 7 (plus a residual error component), but otherwise no 
assumptions of a statistical nature will be required. 


Notation 


It should be noted firstly that the symbols used below to represent 
observations are also used as subscripts to specify the corresponding 
positions in the experimental design. 

Let + represent the set of parameters 7 in the linear model. 

Let z, with elements z, represent a formally complete vector of 
observations for the given experimental design, and let &:2(<), with 
elements &,(«) denote the corresponding vector of expected values 
specified by the linear model. 

Let t denote a set of least squares estimators for the parameters 
7, derived by formally minimizing the error sum of squares 


[z — — 


with respect to t. The ¢ are linear functions of z.’ 

Let Ez , with elements E, , denote the vector of estimators for the 
expected values &,(+). The estimators HZ are obtained by substituting 
the estimators ¢ for the parameters 7 in the expected values &(<). Thus 


E, = &,(t). (1) 


The estimators EZ will be referred to as the plot-estimators for the design: 

In the actual situation to be considered, certain observations are 
missing. Let y denote the vector of existing observations, and let u be 
a vector of unknowns (u, v, w, «+ +) representing the missing values. The 
corresponding vectors of & and E£ will be designated accordingly: &, , &, , 

We have used the term estimator in referring to formulae (functions 
of z). We shall be concerned below with specific estimates, derived by 
applying these estimators to specified sets of data. Thus, for example, 
E,(y, u) denotes the vector of plot-estimates for the missing values, 
derived from the set of data which has been formally completed by 


1Though the point need not concern us here, it should be noted that the set of estimators t may 
not be well-determined. The usual practice is to impose certain linear restraints on the model so that 
well-determined estimators are obtained, but this is largely a matter of convenience. The plot-esti- 
mators EZ, however, are in any case well-determined. 
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inserting the unknowns u, v, w, --- , while £,(y, 0) denotes the vector of 
plot-estimates for u derived from the set of data in which zeros have 
been substituted for the missing values. The latter will be referred 
to as the initial plot-estimates, and the notation will be contracted, 
on occasion, to Ee . 


The missing value procedure 


According to the principle of least squares, the required estimates 
of the parameters 7 from the incomplete data y are those which minimize 
the error sum of squares 


Si(z) = ly — — 

The direct procedure would be to differentiate S,(+) with respect to +, 
but the resulting normal equations would not possess the symmetry 
associated with the complete experimental design, and might in conse- 
quence be difficult to solve. The missing value procedure, on the other 
hand, makes use of the known solutions ¢ = t(z) of the normal equations 
for the complete design. The following argument shows clearly the 
relation between the two procedures: 

Let the data y be formally completed by adding a set of unknowns 
u for the missing values, and consider the minimization, with respect 
to both + and u, of the error sum of squares for the completed data, 


S(2,u) = [y — — + fu — — 

S,() + S.(+, u). 

The minimization of S(z, u) can be carried out stepwise in two ways: 
(a) Minimize with respect to, first u, and then +. Since S,(2) does 

not involve u, S(t, u) is minimized with respect to u for values of u 

satisfying 


ll 


Il 


u= &u(*), (2) 
for then S(*, u), which cannot be negative, vanishes. The partial 
minimum is then S,(), so that this method leads to the direct procedure 
referred to above. 

(b) Minimize with respect to, first +, and then u. As the unknowns 

u formally complete the data, the standard formulae t provide the 
values of < for which S(+, u) is minimized, namely 

t(y, u). (3) 


The partial minimum is, by virtue of (1), 
S(t, u) = [y — Ex(y, u)]’ly — 
+ fu — — w)], 


which is formally the residual sum of squares for the completed data. 


(4) 
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The backward process supplies the missing value procedure suggested by 
Yates [12]: Firstly the data is completed with the estimates u which 
minimize the residual sum of squares (4). The required estimates of 
the parameters 7 are then given by the standard formulae t [equations 


(3)]. 
Equations for the missing values 


It is not actually necessary to differentiate the residual sum of 
squares (4), to obtain equations for the missing values. Equations (2) 
and (3) above, which must be jointly satisfied by the values of + and u 
which minimize S(t, u), provide the necessary equations, for substitution 
from (3) in (2) gives, by virtue of (1), 


u). (5) 

In other words, each unknown should be equated to the corresponding 
plot-estimate, derived from the formally completed data (and thus a 
function of the unknowns). 

The following decomposition shows, in more detail, the structure 
of the equations (5): 

The formally completed vector of observations can be expressed 
as the sum of two components, 


Since the estimators H(z) are linear functions of z, a corresponding 
decomposition of the missing value equations can be obtained, 

Ex(y, 0) + E.(O, u). (6) 
This decomposition separates the numerical from the algebraic terms. 
Further decomposition is effected by expressing the data (0, u) as a 
linear combination of basis vectors 6, , which consist of a unit in the 


z-position and zeros elsewhere. Thus in the case of three missing values 
u, v, and w, 


- 


0 


uU 1 0 
v 0 1 


= ud, vd, + wd. 
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and the corresponding decomposition of the equations is 
u = E,{y,0) + + + wE,(8,) 
v = 0) + + vE,(8,) + wE,(8,) (7) 
w = E,(y, 0) + + vE.(8,) + 


Transferring terms in the unknowns to the left, we may express the 
equations in the form 
Au=E., (8) 
where E: [contracted from E.(y, 0)] denotes the vector of initial plot- 
estimates for the missing values, and A is the matrix of coefficients. 
The diagonal elements of A are a,, = 1 — £,(8,) etc., and the other 
elements are a,, = —£,(8,) etc. Although not obvious in the above 
derivation, it can be otherwise shown that the matrix A is symmetric. 
The computation of the missing value equations will be discussed 
in Section 3. 


Missing values for covariance analysis 


The missing value procedure in this situation has been outlined in 
an earlier paper by the writer [10], and illustrated there with a numerical 
example. Only a formal derivation of the relevant equations for missing 
values will be given here. 

Let z (as before) denote a formally complete vector of observations z, 
and let w denote the corresponding vector of values of a concomitant 
variate w. The appropriaie linear model for this situation specifies 
the expected values 


8) = &.(t) + Bw, (9) 
in which @ is the regression parameter. A formal least squares analysis 


on the basis of this model provides the following estimators for the 
parameters + and for 8, 


t*(z; w) = t(z) — bt(w), (10) 
_ __ Residual 
a Residual §.s. for w’ (11) 


where t, as before, represents the standard estimators for the experi- 
mental design. Substitution from (10) in (9) gives the adjusted plot- 
estimators, 
E%,,,(Z; E,{2) bE.(w) + bw, 
= E,{z) + — E,(w)], 
where FL is the standard plot-estimator. 


(12) 


a 
“oh 
¢ 


er. 


ANALYSIS OF INCOMPLETE DATA 263 

Suppose now that certain observations z are missing. Let y denote 
the vector of existing observations, and x denote the corresponding 
vector of concomitant values w. The values of the concomitant variate 
corresponding originally to the missing observations are irrelevant to the 
analysis of the data y, and may therefore be replaced by arbitrary 
values, determined so as to give the simplest analysis. Accordingly let 
u denote a vector of unknowns for the missing observations z, and let 
v denote a corresponding vector of unknowns for the concomitant 
variate. 

Substituting from (12) in (5), we obtain the following equations 
for u, 

u;Xx, v), 


= E,{y, u) + blv — v)], 


in which b is determined from the formally completed data and is 
therefore a function of both u and v. Clearly these equations will be in 
their simplest form if the arbitrary values v are determined so that 


(18) 


v= E,(x, (14) 
in which case the equations (13) become 


u = E,(y, u). (15) 


In other words, the simplest procedure is fo fit two parallel sets of 
missing values, as was stated in [10]. 

It will be noted that when v is determined in this way, an analysis 
of the sum of squares for the completed concomitant data (x, v) yields 
the correct residual sum of squares for the incomplete data x. Conse- 
quently the variance of the regression coefficient b can be determined 
in the usual way, by the formula 


2 


o 
Vb) = Residual S.s. for 


The analysis of the sum of products for the completed data likewise 
yields the correct residual sum of products for (y, x), since the residuals 
of y and x are correctly given by y — Ey(y, u), x — Ex(x, v). The 
residuals of u and v are of course zero, by virtue of the equations (14) 
and (15). 


Some properties of the matrix of coefficients A 


The plot-estimators H(z) are linear functions of z, the coefficients 
in which involve only the structural parameters of the experimental 
design. (For the reasons given above we may exclude from consideration 
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those estimators E* which give linear adjustment on a concomitant 
variate.) Consequently the matrix of coefficients A in (8) is a function 
of only the structural parameters, and is determined by the positions of 
the missing values in the design. It is apparent, moreover, that the 
value of any particular coefficient a,, , corresponding to a pair of missing 
values (u, v), [the coefficients a,, being understood as corresponding 
to the identical pairs (u, u)], is determined only by the positions of 
u and v in the design, and is not dependent in any way on the positions 
of other missing values. 

The value of a coefficient a,, will depend primarily on the relative 
position, or configuration, of u and v in the design, and only to a secondary 
extent on their actual positions (as in designs with non-uniform replica- 
tion for instance). We shall therefore speak of the coefficient a,, as 
being determined by the relation between u and v in the experimental 
design, the word relation being extended in meaning to include the 
secondary aspect mentioned.” 

The relations between pairs of observations in an experimental 
design are usually of a simple kind, compounded from elementary 
relations such as “‘in the same block” or ‘“‘with the same treatment” 
and the opposite relations ‘‘in different blocks,’ etc. Most experi- 
mental designs in common use possess a high degree of symmetry in 
their structure, in consequence of which the number of possible relations 
between pairs of observations is very limited. It is a simple matter, for 
these designs, to tabulate the various relations and the corresponding 
values of the coefficients a,, which these relations determine. The 
method of constructing such tables, and their application, is explained 
in Section 3. 


3. SETTING UP THE MISSING VALUE EQUATIONS 


Determining the plot-estimators 


The basic formulae required to determine the equations for missing 
values are the plot-estimators for the experimental design, based on 
the appropriate linear model. The required plot-estimators are simply 
the linear combinations, indicated by the model, of the estimators ¢ 
for the parameters r. 

The formulae EF should be multiplied through by a numerical factor 
to eliminate the need for division in the computations, and should be 
expressed in such a form that maximum use will be made of the pre- 
liminary computations which are required, in any case, for the analysis of 


2The idea of relation is only loosely defined here. James [6] has since given the concept more 
precise definition, and has developed the concept of a relationship algebra for an experimental design. 
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variance. In the case of non-orthogonal designs, the estimators for 
certain of the parameters, the estimators for block effects in incomplete 
block designs for instance, may not ordinarily be required for a standard 
analysis. These, however, usually can be expressed simply in terms 
of the estimators normally computed. 

Sometimes, of course, the full linear model, on which the analysis 
of variance is based, is unnecessarily elaborate so far as specifying 
the expected values for particular observations is concerned. For 
example, the decomposition of a total treatment effect into factorial 
components may be irrelevant. It is as well to determine the plot- 
estimators by referring to the simplest equivalent form of the model, 
and thus avoid unnecessary algebra. The most important instance of 
this situation is when the residual sum of squares is a sum of components, 
each derived from a separate region of the experiment. Correspondingly, 
the plot-estimators for any one region are functions only of the obser- 
vations in that region, and will therefore be arrived at most expeditiously 
by considering the ‘local’ model for the particular region. Consider for 
instance the following experiment performed on rats, to investigate the 
effect of two diet factors: 


Treatments: 2 levels of A X 3 levels of B. 

Purpose: To check first for interaction, and then to examine the 
main effects of B. The main effects of A were of no interest. 
Design (split-plot): 10 blocks (whole-plots). Each block consisted of 
4 litters (each of 3 rats) assigned to 3 cages, with four rats, one 
from each litter, in each cage. The treatments A (equally 
replicated) were assigned at random to the blocks, and the 

treatments B, to the cages of each block. 


The skeleton of the analysis of variance is shown in Table 1. In the 
actual experiment the observations on a few rats had to be discarded, for 


TABLE 1 
SKELETON ANALYSIS OF VARIANCE 
(degrees of freedom in parentheses) 
Between blocks (9) A (1) 
Residual (a) (8) 


B (2) 
Between cages (20) {4 xB (2) 
Within blocks (110) | Residual (b) (16) 


Within cages (90) nen (30) 
Residual (c) (60) 
(LE 


| 
|| 
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2 
nds Re 
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accidental reasons. ‘The appropriate procedure was then to estimate 
missing values so as to minimise the residual sum of squares (c). Now 
the latter is the sum of 10 components, each component being the litter 
X cage interaction for a particular block. The appropriate ‘local’ 
model thus specifies merely the additive litter and cage effects for the 
block, and the required estimators are therefore given by the following 
formula, 


12F, = 4L, + 3C. — B., 


L,, C, , and B, denoting the relevant litter and cage totals, and the block 
total, respectively. 
The following examples give further illustration of the procedure 
for determining plot-estimators: 
1. Designs arranged in blocks. The basic formula, assuming additive 
block and treatment effects, is 


E=m+u+4. (16) 
Thus for Randomized Blocks, in the usual notation, 
biE, = bB, + iT, — G. (17) 


In the case of Balanced Incomplete Blocks, block effects are not 
usually determined, but are obtained by adjusting the block mean for 
the treatments which occur in the block: 


k(m + b,) = Dot. 
Bu 
The appropriate formula is therefore 


MkE, = MB, + kQ. — 
By 


Q. = kT. LB. 


The subscript below the summation sign in these formulae indicates the 
field of summation. 


2. Designs arranged with double restriction, in squares or rectangular 
arrays. ‘The usual model gives 


(19) 
Thus for ap X p Latin square, repeated n times, 
np E, = n(pR, + pC. — S.) + (pT. — @), 


(20) 


‘a 
| 
(18) 
where 
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and in particular, ifn = 1, 
pE, = pR, + pC. + pT. — 2G. (20’) 


In a lattice square design, with (k + 1)/2 or fewer replications, the 
row effects are given by the formula 


— 1)r, = —kL, —rS,+G, where L, = 
Ru 


A similar formula, with M instead of L, applies for the column effects. 
The treatment effects are given by 


kr(r — 1I)(m+ + HL+ SM, =Q, say. 
Tw Tw 
(If the design is balanced, Q, simplifies to 
Q=hT,-—r 
Tu 
Combining the various effects, we obtain 
k’r(r — DE, = kQ, — krL, krM, -— + DSL +0 4+ DG. (21) 


Note that when the basic design is repeated, formulae (16) and (19) 
imply additivity of treatment and design-replicate effects. Sometimes, 
however, it will be necessary to allow for interaction between these 
factors, in which case the formula for a single design-replicate should 
be applied to each design-replicate with missing values. 


Methods of computation 
(a) Direct computation of the equations in the formu = E,(y, u). 


This is simply a matter of determining plot-estimates for the missing 
values in a formally complete form, by inserting unknowns for the 
missing values in the data and carrying these forward in the computa- 
tions. Thus each plot-estimate will consist of a numerical term (the 
incomplete plot-estimate) plus a linear combination of the unknowns. 
As a preliminary step, all simple totals of the data should be deter- 
mined. Then, in the case of non-orthogonal designs, various adjusted 
totals will need to be determined. In some cases it will be simpler to 
obtain, in formally complete form, only those adjusted totals that are 
required in computing the plot-estimates. The remainder of the ad- 
justed totals can then be determined in numerically complete form, 
after the missing values have been determined. In other cases it will 
be preferable to determine all adjusted totals in formally complete form, 
in order to use the automatic checks that are usually available. 
Example. ‘The Randomized Block data in Table 2 have three 
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TABLE 2 
RANDOMIZED BLocK EXPERIMENT 


Treatments 
C 


D Total 


34 35 171 
33 v w 34 104+v+w 
35 35 137 + 

30 28 142 


Blocks 


135 


60 +u+ov 


99 + w 554 +u+v+w 


missing values, represented by u, v, and w. The first equation for the 
missing values is 


20u = 4B, + 57, — G, 
= 4(137 + u) + + u +) — (554 +u+v4+ wv), 
= 294 + 8u + 40 — w, 


and similarly for the second and third equations. Transferring terms 
in the unknowns to the left, we finally obtain the equations 


12u — 40+ w = 294 
—4u + 12v — 3w = 162 (22) 
u— 3v+ 12w 


The solution of these equations is discussed in Section 4. 

(b) Using tables of coefficients. If an experimental design is in 
frequent use, it may be desirable to prepare a table of relations for the 
design, and the corresponding formulae for the coefficients in missing 
value equations, as was mentioned in Section 2. If such a table is 
available, the equations can be determined more expeditiously in the 
form 


Au = Ex, 


the matrix A being determined by reference to the table of coefficients. 
The initial plot-estimates are computed simply from the incomplete 
totals of the data. 

The method of constructing such tables will be illustrated for the 
simple Randomized Blocks design. Consider a basis vector 5, , repre- 
senting a set of observations which are all zero except for a unit in the 
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plot designated v. Let u designate any plot in the design. The block 
and treatment totals relating to the plot wu, and the grand total, derived 
from 6, , are as follows: 


B, 


1 if w is in the same block as », 


0 otherwise. 


= 1 if u has the same treatment as », 


0 otherwise. 
G, = 1. 


Substituting in the plot-estimator for u, given by 


bik, = bB, + tT, — G, 
we find that 


bt =b+t-—1 if u =», 
=b-1 if u is in the same block as 2, 
={— 1 if u has the same treatment as », 
= otherwise. 


The corresponding table of coefficients for the design is therefore 
[referring to the definitions given with (8)], 


bt du» 


(u, v) relation | Same treatment | Different treatments 


1) 


Same block | (b —1)(t — 1) 
1 


Diff. blocks —(t — 1) 


Tables for several of the standard designs have been prepared by the 
writer, and are presented, in compact form, in Tables 3 and 4. It should 
be noted that, in these tables, the part common to the several relations 
in any given line is shown in the marginal column, and the remaining 
part of each relation is given in the body of the table, together with the 
corresponding coefficient. Codes are used to specify some of the more 
complex relations, the explanations of which are given in the blank 
areas of the tables. 

The method of constructing a matrix of coefficients A from these 
tables is straightforward: The relation between each pair of missing 
values is determined by inspection of the experimental design. The 
corresponding coefficient in the table is then evaluated and inserted 
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in the matrix. (An alternative procedure is to evaluate all the coeffi- 
cients in the particular table, before proceeding to the construction 
of A.) 

Examples. Consider the configuration of missing values in Table 1. 
The relations between the missing values are as follows: 


u v w 

u Identical Same T Diff. B and T 
v Identical Same B 

w Identical 


The corresponding matrix of coefficients, from the table for Randomized 
Blocks, is therefore 


12 —4 1 
1-38 


as obtained previously in (22). 

To take a more complex example, consider the data given in Table 5 
(data from Cornish [2]). The design is a balanced incomplete block 
design, with b = 12,t = 9,k = 3,r = 4,\ = 1,N = 36. The treat- 
ments are lettered a to 7, and five observations assumed missing are 
designated u, v, w, xz, and y. It can be seen firstly, that no two missing 
values have the same treatment. Secondly, only the pair (u, v) occurs 
in the same block. For the other pairs of missing values, the numbers 
(s) of treatments common to the corresponding pairs of blocks, are 
given below. It is also necessary to note in each case the value of q 
(defined in Table 3). The non-zero values of q are shown in parentheses: 


(s) v w y 

u (Same B) 0 0 1 

v 0 0 1 

w 0 1(q = 1) 
x 1(q = 1) 


With this information the matrix A can now be determined by reference 
to Table 3. The initial plot-estimates E° are given by the formulae 


27E, = 9B, + — Q, ete., 


By 


\ 
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TABLE 5 


BALANCED INCOMPLETE BLocK EXPERIMENT 


a u d w g = b  =66.0 a 61.0 c 69.9 
b «666.2 e 63.9 h 63.1 ad f d 62.0 
c v f 67.8 ¢ 65.3 t y h 62.9 h_ 62.6 
66.2 131.7 128.4 127.1 190.9 194.5 
+u+v +w +z +y 
65.4 c 69.2 b 66.6 a 56.9 58.8 
e 64.8 e 64.3 d 61.9 63.4 
h_ 638.9 g 60.2 g 61.0 g 61.9 t 65.8 & 65.1 
194.1 193.7 194.8 180.7 202.2 187.3 
a b c d e i g h t Total 


T| 176.7 264.2 208.2 185.0 256.4269.3 183.1 2525 196.2 | 1991.6 
+u +0 +w +2 +y +u+... 
Q| —95.0 210.4 -—32.0 —-—79.0 62.4 88.3 —148.3 49.6 —56.4 0 
+2u  —u-v—y —ut+2v 4+2w-y —w +22 +0 


For the block containing: 
uandv w y 


@ m4 V4.9 —155.1 75.0 


and the equations finally obtained (multiplied by the factor 27) are 

as follows: 

[12 -6 0 O [227.4] 

12 0 » 416.4 
0 O12 =| 876.6 }- (23) 
0 0 0 12 22iiz 865.8 

—-1 2.2 | 899.7 _| 


The solution of these equations will be discussed in Section 4. 
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4. SOLUTION BY MATRIX INVERSION 
The equations for missing values are usually, though not always, 
non-singular. We shall consider first the non-singular case, and deal 
later with the singular case. 
Some inversion formulae 


The solutions of the equations Au = Zi may be computed in the 
form 


u=-A'E, 


in which A“ is the inverse of the matrix A. The most useful formulae 
for the inversion of non-singular symmetric matrices are listed below: 


(a) A |A| [A;;], (24) 


in which the A,; are the cofactors of the elements a,; of A. 
(b) [A ia + PQP’ 
Cc’ B — QP’ Q 
where P = A''C, Q = (B — C’P)”. 
u(u + 
where I is the unit matrix, and 1 denotes a column of units. 


1+ AD 


(d) (B+ A11’)" = B" 33’, (if B non-singular), (27) 
where 6 is the column vector of row sums of the elements in B™’, and 
D is the sum of all elements in B™*. 


(It should be mentioned that slight modification of (24), (25), and (27) 
is necessary to render these formulae applicable to non-symmetric 
matrices. ) 

Since the matrix of coefficients A, when multiplied by the appropriate 
scale factor, consists only of simple integers, the application of these 
formulae will require, to a large extent, only mental arithmetic. 


Applications 


(a) Two or three missing values. For two missing values we have 


c b —¢ 


ate 
: 
¢ 
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For the case of three missing values, formula (24) should be used unless 
(26) is applicable. Each cofactor A;; is determined by deleting (men- 
tally) the row and column in A containing a;; , and computing the cross- 
product for the remaining 2 X 2 matrix, multiplied by the sign factor 
(—)‘*’. The determinant |A| may be determined by computing the 
sum of products of any row of A with the corresponding row of [A,,]. 
For example, referring back to the equations (8), we have 


12 —4 1 


135 45 
[Ai] =| 45 143 32], 
0 32 128 


|A| = 12 X 185 — 4 X 45 = 1440. 
The solutions of the equations are therefore 


u = (135 X 294 + 45 X 162 + 0 X 357)/1440 = 32.625 


and, similarly, v = 32.208, w = 35.333. 


(b) Four or more missing values. Generally the formula (27), and 
sometimes (26), will be applicable (see below), but otherwise it will 
usually be best to partition the matrix suitably, and apply formula 
(25), as in this way the simple, integral character of the matrix can be 
largely retained in the computations. 

Consider, for example, the equations (23). Here the submatrix in 
the first four rows and columns can be inverted at sight, and in the 
notation of formula (25), 


0 


0 
2 
0 


B-CP=l11, Q=1/ll, 


| 

4 12 -3 | 
1 -3 

“i, 
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| 

0 0 anal 2400 

36 0 0 | ‘ 

0 3 I 
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1 -1 
| — 


1 
-PQ=— 


The full inverse is therefore 
45 
23 
| 

6 


= 


which gives the solutions 

u = 59.285, v = 69.758, w = 62.192, x = 61.192, y = 65.145. 

If the inverse matrix is not otherwise required, the solutions of the 
equations in the partitioned form 


can be determined directly by the formulae 
v = QE, — 
u = — Pv. 


(28) 


(c) Missing values with a single relation. With only a few missing 
values, it will often happen that they have the same relation to each 
other. They may be all in the same block for instance. More usually 
they are isolated, each in a different region or replicate of the design. 
In these cases formula (26) is applicable, and the solutions of the equa- 
tions Au = Eu (for r missing values say) are given by the formulae 


0 
LE » etc. (29) 

(d) Grouped missing values. If there is more than one relation 
between the missing values, it is frequently possible to classify them into 
groups so that the missing values in any one group have the same relation 
to the missing values in any other group. Correspondingly, the matrix 
of coefficients can be expressed in the form A = B + A1I1’, in which B 
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23 -1 -1 6 
45 -1 -1 6 
—1 34 1 
—1 1 34 
6 -6 -6 36 
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is the direct sum of submatrices B; , each corresponding to a group of 
missing values. In other words, with an appropriate ordering of the 
missing values, 


B, 0 0O 
0 O 
oO 


B= 


Application of the formula (27) reduces the inversion of A to the inver- 
sion of the matrices B,; (which will usually be only of order 1, 2, or 3), 
and provides the following solutions of the missing value equations, 


i=1,2,---, 


(30) 


where u, (sub-vector of u) is the vector of missing values in the 7-th 
group, 5; (sub-vector of 6) is the vector of row-sums of B;', and W = 
s’Et . If w, and d,; denote the quantities for the 7-th group, corre- 
sponding to W and D, then W = ,D = d;. 

An important application of the foregoing is to the estimation of 
missing values in Randomized Blocks, which will be considered below. 

It will sometimes happen that an otherwise appropriate grouping 
of the missing values is spoilt by one or two links (usually by treatment) 
between missing values in different groups. In this case the trouble- 
some elements in the matrix should be partitioned off, so that the above 
considerations apply to the remaining submatrix [when applying formula 


(25)). 


Estimation of missing values in Randomized Blocks 


We shall avoid using the specific terms blocks and treatments below, 
and shall refer instead to the p rows and g columns of a double classifica- 
tion. 

The appropriate grouping, for the application of formula (27), of 
the missing vlues in this p X q classification, is into groups which 
have no rows or columns in common with each other, and the intergroup 
relation is then ‘in different rows and columns.’ Correspondingly, if 
A denotes the matrix of coefficients multiplied by the factor pq, 


A=B+11’, 


where B is the direct sum of matrices B; , the elements of which are 
given by the following table: 
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Dur 
(u, v) relation Same C Diff. C 
Same Rk 
Diff. R —q 0 


Let n denote the vector pgiu . The elements of n are then 
m = pR, + pC. — G, etc., 


where R, , C, , and G are the relevant row and column totals and the 
grand total of the incomplete data. In this notation the missing value 
equations are Au = n, and their solutions are given by the formulae 
similar to (30), (with A = 1), 


W 
u; = B;'n; — 
1+D (31) 
W 


in which the quantities u* will be referred to as the incomplete estimates. 

For a given group of missing values, the incomplete estimates u* 
and the associated quanties 5; , d; , and w,; , will depend on the con- 
figuration of these missing values in the p X gq classification. For the 
most commonly occurring configurations, the necessary formulae for 
these quantities are given in Table 6. The table also gives, with each 
configuration, the complete estimates u for the special case in which 
all missing values are in the one configuration; that is, when there is 
only the one group. Configuration (1) in Table 6, incidentally, really 
corresponds to / groups, each with only one missing value, but for 
practical purposes these isolated values are best considered as one 
group. 

Table 6 covers some special cases considered earlier by Thompson 
[8]. It should be noted that the formulae in Table 6 contain implicitly, 
and in a concise form, the necessary information for determining the 
inverse of A, if this is required. 

To illustrate the application of Table 6, we shall consider the example 
given in Table 2 of Yates [12]. There are 8 treatments, 10 blocks, and 
9 missing values which, with some rearrangement of blocks and treat- 
ments, have the following configuration: 


pace 
mele 
| 
| 
Ny 
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Blocks 
136578 
NK |w | 
Treatments 0 |- 
NKP ue! 
KP 
NP 
N its * | 


Our notation for the missing values differs from Yates’ notation, the 
correspondence being as follows: 


U2 Us Us Us Us Uz Us Ug 


a b e c d g a f h 


The configuration shown above breaks up into three disjoint configura- 
tions, corresponding to 1, 4, and 5 transposed, in Table 6. 

The first stage of computation is to determine the column of quan- 
tities ». The values given in Table 7 have been computed from the 
incomplete totals of Yates’ data (not given here). The preliminary 
quantities required in Table 6 are then determined: 


1. p = 8,p’ = 7,¢ = 10, q’ = 9. a = 62, S = 399.14. 
4. p=6,0=8. q'p = 54, p’o = 56, po = 48. 
a = 432,868 = 560, 7 = 2976. = 62,e = 4960. 
l = 172, a = 2852. P = 431.74, Q = 399.42. 
5. (w.r.t. transposed classification.) p= 10, p’ = 9, p” = 8,q = 8. 
1=n=1,0, =a =p; =p =}. 
a= 106. m=17,b= 203. Q 357.12,» = 0.918. 


The quantities w are also determined at this stage, and are recorded 
under the sub-columns of n in Table 7. The next column 6 in Table 7, 
with sub-totals d, is then entered in quotient form. The total of the 
column 6 is computed as D = 2/62 + 172/2852 + 17/203 = 0.1763. 
Finally the factor W/(i + D) is computed: W/(1 + D) = (6.438 + 
12.650 + 16.563)/1.1763 = 30.31. 
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TABLE 7 
CoMPUTATION OF MissING VALUES 


No. Configuration n 7 u 
1 uw 209.11 1 2.884 
+ 62 
2 190.03 1 2.576 


w = 6.438 d =2 + 62 


3 200 .07 62 3.757 
4 231.67 54 + 2852 3.733 


5 199.35 56 3.333 


w= 12.650 | d = 172 + 2852 


6 195.01 4.5 3.606 

7 us ug F 162.11 4.5 3.218 
+ 203 

8 Uz ug 199.73 4 3.314 

9 239 .07 4 3.886 


w = 16.563 d = 17 + 203 


D = 0.1763 Wil + D) = 30.31 


The estimates u of missing values may now be computed in one step. 
In each case, first compute the incomplete estimate u* (except for u, and 
Uz), using the appropriate formula from Table 6, transfer to the key- 
board and multiply by the divisor in the corresponding 6, substract the 
product of W/(1 + D) with the numerator of 6, and divide by the 
divisor of 6. Typical calculations are as follows: 


u, = (209.11 — 30.31)/62 = 2.884. 

uf = (8 X 231.67 + 431.74 + 10 X 12.65)/560, 

U, = (uf X 2852 — 30.31 X 62)/2852 = 3.733. 
uz = (199.73 — 0.918 + 5 X 16.563)/72, 


Us = (u¥ X 203 — 30.31 X 4)/203 = 3.314. 


- The final estimates should be checked by re-computing the quan- 
tities n using completed totals, the check relation then being u = »,/pg. 
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Comparison with iterative method. A count of the operations of 
setting up, adding, multiplying, recording, etc., and of the digits in- 
volved, shows that in the above example, the computations are equiv- 
alent to two cycles of iteration with the formula for a single missing 
value. Checking is equivalent to an additional cycle. 

The example, however, is a most unfavorable one for comparison, 
with seven out of nine missing values in the two most complicated 
configurations tabulated. With more favourable configuration, the 
amount of computation would be more of the order of one and a half 
cycles of iteration, and in the simplest cases, in which all missing values 
are in one configuration of the type 1 or 2, the effort involved is little 
more than for a single cycle. Moreover, the computations will give 
any required degree of accuracy. 

Against this, the iterative method usually requires two, sometimes 
three, or even four cycles to obtain sufficient accuracy, plus an additional 
check cycle. 


Solution of singular equations for missing values 


Singular equations Au = n, if consistent, have an infinite number 
of solutions, each of which can be expressed in the form u = Cn, where 
C is a matrix satisfying the relation 


ACA = A; (32) 


for if u, is some solution of the equations, we have ACn = ACAu, = 
Au, = n (James [5]). The matrices C, of which there are an infinite 
number, are the effective inverses of A.* 

Singular equations for missing values arise when whole sections, 
such as blocks, of a design are missing. It is clear from the derivation 
in Section 2, that any set of estimates satisfying the missing value 
equations will lead to a correct analysis of the data, so that the choice 
is only a matter of convenience. 

In practice, inspection of the matrix A will indicate a symmetric, 
idempotent matrix F (that is, a matrix F such that F’ = F), and a 
non-singular symmetric matrix B, such that 


A = BF. (33) 


The matrix F is uniquely determined by A, and is, in fact, the matrix 
which gives the orthogonal projection, of any vector u, on the range of 
A (that is, on the space of vectors Au). The matrix B, however, will 
be to some extent arbitrary. The inverse B~ of B is clearly an effective 


3A generalized inverse for matrices was defined by Penrose [7]. This inverse is unique, and deter- 
mined by four relations. An effective inverse (as defined above) satisfies one of these relations. 


| 
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inverse of A, and gives the solution which lies in the range of A, 
uo = Fu = B'n. (34) 


Any other solution u of the equations can be expressed in the form 
u = uy + Vv, where vis any arbitrary vector orthogonal to uy , or that is, 
such that Fv = 0. Thus Fu represents the estimable component of u, 
and the linear functions orthogonal to Fu, represented by Gu, where 
G = I — F and hence FG = 0, may be arbitrarily determined. 
Examples. Consider a B.1.B. design with one missing block (k plots). 
The equations for the missing values may be simplified to the form 


(At — — 


or in vector form 


Il 


Ni» +=1,2,---k, 


k)Fu=n 
where F = I — 11’/k. (F’ = F). The estimable component of u 
therefore has elements 
Ui =U; — = /(At — += 1,2, --- k. 
G = 11’/k determines the linear function a, which may be given an 
arbitrary value uw. The general solution is therefore 
u; = n./(At — += 1,2,---k. 


To take a more complex example, suppose that, in the B.I.B. data 
of Table 5, all of the first block is missing, the three missing values being 
u, v, and v’ say. The matrix of coefficients for the six missing values 
(five previously) can be expressed in the form (33) as follows: 


A B F 

12-6 -6 0 O © 0 =1 0 0 O 
-6 12-6 0 0 2 0 2 3-3 0 0 
-6 12 0 O -1]= 0 0 18 0 O -1 #0 0 0 
0 0 012 0 2 0 0 01 0 0 
-1 2-1 2 2 12 -1 2-1 2 2 12. 0 0 00 0 1 


The matrix B can be inverted easily, using formula (25), and hence the 
estimates uy determined. Alternatively the equations for uy = Fu may 
be solved directly by the formulae (28). As in the previous example, 
the block mean for the missing block may be arbitrarily determined. 
The missing values w, x, and y, however, will be well-determined. 
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QUERIES AND NOTES 


GrorcE W. SNEDECcOR, Editor 


130 QUERY: Orthogonal Coefficients for Unequal Intervals 


In an experiment performed here, I used four levels of a nutrient: 
0, 5, 10,20 mg. I cannot find a method for calculating the three sets of 
orthogonal coefficients for these unequal intervals. If it can be done, 
please advise me. 


A discussion of this matter, together with references to 
ANSWER: other work, was given by Wishart and Metakides, ‘“‘Orthog- 

onal Polynomial Fitting,” in Biometrika 40, pp. 361-369 
[1953]. The method given allows not only unequal intervals but also 
different weightings at the several levels. 

Some workers may not be familiar with the abbreviated schemes for 
computing matrices. If the number of levels is small, as in your experi- 
ment, or if regressions up to only the 3rd degree are needed, a less 
formal method is rather easily applied. It will be illustrated by use of 
your levels, coded to 0, 1, 2, 4. 

In Table I the linear coefficients are symbolized by §,; = X; + a. 


TABLE I 
CoMPUTATION OF COEFFICIENTS FOR LINEAR REGRESSION 
X; =X: +a ts 
(1) (2) (3) (4) 
0 a —7/4 -—7 
1 1 + ay —3/4 
2 2+a4 1/4 1 
4 4+a,. 9/4 9 
Sum 7+ 4a, =0 0 0 
= —7/4 
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The values of X; in column (1) are substituted successively in £,,; , the 
results being recorded in column (2). One of the fundamental properties 
of the £; is that their sum is zero. This leads toa, = —7/4. Substitu- 
tion of this value of a, in column (2) gives the £,; of column (3). To 
reduce those to the simplest integers, multiply by 4; the resulting ¢/ 
coefficients are in the last column. In many cases, these will be the 
only coefficients required. 

If coefficients ~; are needed for the second degree regression, set 
them equal to Xj + b.X,; + a, , then substitute the values of X; . 
The results are in Table II, column (3). The sum of these coefficients 


TABLE II 
CoMPUTATION OF COEFFICIENTS FOR QUADRATIC REGRESSION 
Xi thi for = WX: + a2 (Xi = bX; + a2) £35 £4; 
(1) | (2) (3) (4) (5) (6) 
0 -—7 a2 — 2 7 
1 1 + be + a: -3- 3d2 —8/7 
2 1 4 + 2b. + ae 4 + 2be + a2 —16/7 -8 
4 9 16 + 4b2 + ae | 144 + 36b2 + 9a2 10/7 5 
Sum 0 21 + 7b. + 4a, = 0 145 + 35b = 0 0 0 
—145/35 = —29/7 


~ 
w 


—1/4 {21 + 7(-—29/7)} = 2 


must be zero, as before. Now apply the second fundamental property 
of the coefficients; the sum of the products, &,;&,; (or &{,&;) must be 
zero. These products, column (2) by column (3), are recorded in 
column (4); their sum is set equal to zero and solved for b, (= —29/7). 
This result, substituted in the sum of column (3), leads to a, = 2. 
a, and b, are now substituted in the £.; of column (3) with the results 
shown in column (5). Multiplying by 7/2 yields the ££, of the last 
column. 

In practice, the coefficients for the third degree regression are rarely 
needed; the sums of squares for the first and second degree regressions 
are subtracted from the total sum of squares to get the remaining sum 
of squares for the third degree regression. If the coefficients are required, 
the formula, §&; = X} + c,X{ + bX; + a; , is manipulated in the 
manner of Table II. There are now three relationships: > &; = 0, 


wt 
we 
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> fits = 0, > &,é. = 0. The three equations are: 


73 + 2ic; + 7b, + 4a; = 0 
581 4- 145c, + 35b; = 0 


252 + 44c; = 0. 


Solving: c) = —63/11, b; = 392/55, a, = —36/55. 
Substituting these values in the &;; : 


= — 36/55, = 96/55, 
Multiplying by 55/12, 


= —3, = 8, = —6, 


— 72/55, = 12/55. 


=1. 


It is clear that with more levels the calculation of all the coefficients 
will become tedious. Those familiar with the methods used by Wishart 
and Metakides will prefer that scheme of computation. 
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ABSTRACTS 


The following papers were presented at meetings of the Biometric Society 
and the American Association for the Advancement of Science in India- 
napolis, Dec. 27 and 28, 1957. 


JOHN 8. CHIPMAN AND M. M. RAO (University of Minnesota, 
487 Minneapolis, Minnesota). On the Use of Idempotent Matrices in 
the Treatment of Linear Restrictions in Regression Analysis. 


It is desired (I) to estimate the kth-order column vector 8 in the 
regression equation 


y= (1) 


where y and e are nth-order column vectors of n observations on the 
dependent variable and true residual, and X is the n X k matrix of 
observations on k independent variables, n > k, subject to the restriction 


WB =a (2) 


where WY is a fixed r X k matrix of rank r < k and a a fixed column 
vector of order r. (II) The hypothesis (2) is tested given (1) and the 
assumption that X is fixed in repeated samples and ¢ has the multi- 
variate normal distribution with E(e) = 0 and variance matrix 
V(e) = Io’. Finally (III) the general linear hypothesis is considered, 
in the form of testing one set of linear restrictions V,8 = a, subject to 
another set "8 = a, being true, where VY , VW, are ro X n of rank ro 
and r,; X n of rank r, respectively. | 

The transformation 78 = §* is defined, with T’ = [U’W’], where 
U is some 1 X k matrix (1 = k — r) such that T is nonsingular, and ’ 
denotes transposition. The inverse transformation is denoted by 
T~' = [®'V"], a partition into 1 and r columns. The matrices L = &’U 
and R = V’W are idempotent of ranks / and r respectively, and orthog- 
onal complements (LL = L, RR = R, LR = RL = 0). Thus they 
effect projections of the k-dimensional parameter space on a linear 
manifold £ along a linear manifold ®, and vice versa, and the effect of 
the restriction (2) is to fix an element p « ® and confine 8 to the hyper- 
plane £ + {p}. 
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The solution to (I) is given by 
B = b — "(Wb — a) (3) 


where b = (X’X)~'X’y is the unrestricted least-squares estimate of f. 
The restricted estimate f of 8 is given by a projection of the unrestricted 
estimate b of 6 on the hyperplane £ + {:V’a} along a linear manifold @, 

A test (II) of the hypothesis (2) is given by the variance ratio 


n—kéHe — — a). 
r eFe yy — y Xb 


(4) 


which has the F distribution with r and n — k degrees of freedom, 
where H = X(X'X)'W[W(X'X) W(X'X) 7X’ and E = — 
X(X'X)"'X’ are symmetric idempotent matrices of ranks r and n — k 
respectively, and HE = EH = 0. Thus H and E£ effect perpendicular 
projections of the n-dimensional sample space on linear manifolds of 
dimension r and n — k respectively, and Cochran’s theorem follows 
from the Pythagorean property @’é = e’e + h’h, wheree = Ee,h = He, 
and y — XB = (E+ Ade. 

A natural extension of the preceding analysis leads to a solution 
of (IID), given by the F ratio 


n—-k+7 
+ (5) 
where Hy = *X’, H, = H— Hy, 
and W’ = [W/W]. 


DONALD R. HILL (Minnesota Mining & Manufacturing Com- 
488 pany, St. Paul, Minnesota). Some Uses of Statistical Analysis in 
Classifying Races of the American Shad (Alosa Sapidissima). 


Each year pound nets fished in the ocean off the coast of New York 
and New Jersey catch large quantities of shad. Most of these fish are 
native to the Hudson and Connecticut Rivers and should be considered 
in any management plan for the two rivers. However, estimates of 
the racial composition of this catch are required before these catches 
can be included. 

An analysis of variance of some meristic counts made on these fish 
indicated that the differences between rivers were considerably larger 
than the differences between years within a given river. This informa- 
tion provided additional evidence for the existence of races. 

Using these meristic counts, a discriminant function was constructed 
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which would separate a large proportion of a mixed population. The 
function was constructed from data obtained on Hudson River shad 
in 1939 and Connecticut River shad in 1945. 

The meristic data from a sample of 105 fish obtained from the 
Hudson River in 1940 were substituted into the discriminant function. 
Very good agreement was obtained between the theoretical number of 
misclassifications and the actual number. Methods of estimating the 
racial composition using the classifications obtained from the dis- 
criminant function are also discussed. 


EUGENE A. JOHNSON (University of Minnesota, Minneapolis, 
489 Minnesota). Monitoring and Evaluating Treatment Effects in 
Epileptics by a Graphical Sequential Test. 


Seizures in certain epileptics may be recorded cumulatively with 
time on a graph delineated by two boundary lines. These lines are 
computed (following Dvoretsky, Ann. of Math. Stat. 24, p. 254, 1953) 
with the assumption of a Poisson distribution for time patterns for 
seizures. If the cumulative seizure plot, made by non-professional 
personnel, crosses one of the boundaries the physician is notified. The 
procedure thus aids the physician as a monitoring system yielding an 
analytical statistical appraisal of the relation of current to past seizure 
incidence (i.e., a decision as to the value of a change in medication in 
clinical practice or research). Boundary lines so set as to permit 
decisions “no change or worse” and “improvement” are advocated for 
most epileptics. The simplicity of computation and plotting is empha- 
sized. The empirical basis for applying the procedure + seizure data 
is discussed and the yalue of the procedure for selected cases is illustrated. 


D. KODLIN (University of Pittsburgh, Pittsburgh, Pennsylvania). 
490 Estimation of Risk When Units are Sacrificed Periodically During 
the Follow-up Interval. 


In certain types of experiments the treatment response cannot be 
observed unless the unit is “‘sacrificed.’”’ For instance, mice treated 
with substances producing lung tumors must be killed for diagnosis. 
Though the diagnostician is unable to determine the actual appearance 
time of the tumor, if he finds one, he wishes to estimate the parameters 
of the appearance time distribution after a certain follow-up time during 
which he has sacrificed and diagnosed subsamples. Also required is 
an estimate of the proportion ‘‘susceptible,” i.e., the proportion eventu- 
ally responding. 

The appearance time distribution is assumed to be log normal and 
maximum likelihood estimates of the three parameters are obtained by 
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iteration procedures, data by Andervont, Grady, and Edwards being 
used for illustration. 


BASILIO ROJAS (Iowa State College, Ames, Iowa). The Analysis 


ate of Groups of Similar Experiments. 


A problem of common occurrence in agronomic research is the 
interpretation of a group of experiments consisting of replicated trials 
repeated over locations and years. Similar problems occur in other 
research activities. 

A mathematical model dealing with the analysis of a group of experi- 
ments is proposed which includes many of the models for such analysis 
in the statistical literature as special cases. This model shows that the 
tests of significance offered by other models are unbiased only under 
some assumptions which are unduly restrictive under certain experi- 
mental situations. Tests of significance and estimation procedures are 
given for different models involving various combinations of random and 
fixed components. The common case of combining experiments with 
unequal numbers of replicates is taken up and the bias of the tests of 
significance in the unweighted means analysis is examined. 


R. W. TOUCHBERRY (University of Illinois, Urbana, Illinois). 
492 The Use of an Inverse Matrix in the Analysis of Crossbreeding 
Data. 


Much of the data resulting from the Illinois project on crossbreeding 

dairy cattle is described by the model 
= m+ a; + +e, + acy, + + 

in which m is the mean, a; the effect of the 7-th breed to the k-th breed of 
dam, and eé;;,; is the error variation associated with the /-th calf of the 
ijk-th breeding group. Various characteristics such as weight, height 
at withers, and circumference of chest have been recorded on all project 
animals at six-month intervals from 6 to 48 months of age. The nature 
of these data makes it possible to use one inverse matrix in analyzing 
the data on all characteristics at a given age and in some cases the same 
inverse matrix will suffice for all characteristics at a number of ages. 
The results of the analysis will be given in the paper. 


ROBERT F. WHITE (Iowa State College, Ames, Iowa). Cumu- 


= latively Grouped Response Times in Quantal Response Data. 


A procedure is given for analyzing quantal response data where 
groups of individuals have been exposed to doses of a stimulus and 
observed at pre-chosen times. For each individual, it has been deter- 
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mined that it responded between two observation times or that it 
survived the last observation time. This is more convenient than 
continuous observation to determine exact response times, and more 
informative than observing each individual only once, so that if it 
responded, only the upper limit of its response time is known. In this 
latter procedure the observed responses at the various dose-time com- 
binations are independent, but in the present procedure the observed 
response at each dose-time combination is an accumulation of the 
observed responses at the preceding times. The observations are 
therefore not independent within doses and this increases the complexity 
of the analysis, especially if maximum likelihood estimation is used. 
However, minimum modified chi-square estimation procedures can 
handle the analysis without great difficulty. The aim is to express 
proportional response as a function of dose and time, so that the dose 
which causes a response of r% in given time ¢[ED,(t)] or the time which 
causes a response of r% with given dose d[ET,(d)] can be estimated 
for any d, t. 
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THE BIOMETRIC SOCIETY 


International 


C. H. Goulden of Ottawa, Canada, has been elected President of 
the Society for 1958. The new Council members are: W. U. Behrens, 
A. Buzzati-Traverso, G. Darmois, W. J. Dixon, K. Mather, C. R. Rao, 
and J. W. Tukey. 


Membership 


The paid-up membership of the Society at the end of 1957 was 
1410, an increase of 8% over 1956. Regional figures were: 


Total Total 
Australasian 69 India 14 
Belgian 59 Italian 73 
Brazilian 70 Japan 45 
British 181 Netherlands 42 
Denmark 14 Sweden 14 
ENAR 479 Switzerland 22 
French : 64 WNAR 102 
German 97 At Large 51 


International Standards Organization 


E. van der Laan represented the Society at the third meeting of 
Technical Committee 69 on Statistical Treatment of Series of Observa- 
tions held in the Hague, 15-18 October, 1957. Several other members of 
the Society were included in national delegations. 


Région Belgique et Congo Belge 


The following papers have been read at recent meetings: Factorial 
methods in the analysis of experimental results (Martin); Response— 
surface analysis of a factorial experiment on sugar-beet (Lenger); 
Statistics in the pharmaceutical sciences (Bontemps); Assay of amylo- 
caine (Deltombe). The first two papers were read at joint meetings 
with the Association des Ingénieurs sortis de l'Institut Agronomique 
de l’Etat, the others at a joint meeting with the Soc. Belge des Sciences 
Pharmaceutiques. 
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The following hold office for 1958: President: R. Laurent; Secretary: 
L. Martin; Treasurer: A. Rotti. 


British Region 


The following papers have been read at recent meetings: The struggle 
for existence, with applications to Tribolum populations (E. L. Scott, 
WNAR); Relative growth of interacting organ systems (Skellam); 
A new cross-over design (McDonald); Dynamic of fish populations 
(Pope); Equilibria in Autotetraploids (Parsons). 

The Region collaborated with four other societies in a very suc- 
cessful two-day symposium on Quantitative Methods in Human 
Pharmacology and Therapeutics. J. H. Gaddum introduced the Sym- 
posium, and papers read by members included: Newer statistical 
methods (W. G. Cochran, ENAR); Incomplete blocks in an oxytocic 
assay (Schild); 9 X 9 Latin square in an anti-histamine assay (Dare); 
Analysis of semi-quantitative data (Spicer); Small scale drug trials 
(Bain). 

The officers for 1958 are: President: J. O. Irwin; Secretary: C. C. 
Spicer; Treasurer: P. A. Young. 


ENAR 


The Spring Meetings of the Biometric Society and the Institute 
of Mathematical Statistics were held at Gatlinburg, Tennessee, under 
the auspices of the Oak Ridge National Laboratory, April 10, 11, and 
12, 1958. Thirty-two papers were presented at these meetings. Ab- 
stracts will appear in the coming issue of Biometrics. The following 
attended: G. E. Albert, Mrs. G. E. Albert, G. N. Alexander, Willard 
O. Ash, G. J. Atta, Earl Atwood, R. E. Bargmann, V. P. Bhapkar, 
Allan Birnbaum, Barbara Bishop, C. I. Bliss, V. J. Bofinger, R. C. 
Bose, G. E. P. Box, R. A. Bradley, A. E. Brandt, L. 8. Brenna, W. N. 
Carey, Jr., Charles Carroll, D. Chaudhuri, Victor Chew, Mary Ann 
Cipolloni, A. Bruce Clarke, W. H. Clatworthy, C. W. Clunies-Ross, 
A. C. Cohen, W. D. Commins, W. E. Conner, William 8. Connor, 
Dennis Cooke, Richard Cornell, Jerome Cornfield, L. C. A. Corsten, 
Constance Cox, Edwin L. Cox, Gertrude M. Cox, J. B. Cox, Elliot M. 
Cramer, Jonas M. Dalton, H. A. David, Earl Diamond, N. Draper, J. R. 
Duffett, David Duncan, Arthur M. Dutton, Lila Elveback, Jean Engler, 
A. L. Finkner, Jack Fleischer, $8. M. Free, R. J. Freund, D. A. Gardiner, 
John J. Gart, A. Garzadela, E. Gehan, 8. Geisser, H. Ginsburg, W. 
Glenn, R. Gnanadesikian, Arnold Grandage, J. E. Grizzle, Joan M. 
Gurian, M. A. Guzman, R. J. Hader, W. Haenszel, W. J. Hall, Eugene 
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K. Harris, Boyd Harshbarger, Hubert Hill, W. H. Horton, Wassily 
Hoeffding, J. F. Hudson, R. G. Hoffman, David Hurst, Harold Hotelling, 
Paul E. Irick, M. A. Kastenbaum, Therese Kelleher, George Kennedy, 
Allyn W. Kimball, Marcus Kjelsberg, Carl F. Kossack, Roy R. Kuebler, 
R. G. Laha, E. L. LeClerg, G. J. Levenbach, H. L. Lucas, Eugene 
Lukacs, Mary Lum, Nathan Mantel, Frank Martin, Margaret P. 
Martin, P. A. Miller, D. F. Morrison, George Morton, V. K. Murthy, 
M. D. Nefzger, G. E. Nicholson, Paul S. Olmstead, D. Quade, H. F. 
Robinson, A. C. Rohloff, J. B. Roy, Chas. F. Sarle, M. A. Schneiderman, 
Oliver A. Shaw, 8. S. Shrikhande, Paul Somerville, D. E. South, Harold 
Storz, R. J. Taylor, G. W. Thomson, Malcolm Turner, R. E. Walpole, 
G. 8. Watson, M. B. Wilk, E. J. Williams, R. L. Wine. 


Région Frangaise 


La Société Francaise de Biométrie, 4 la derniére Assemblée Générale, 
a procédé aux élections suivantes: Président: M. Vessereau; Secrétaire- 
Trésorier: M. Schwartz (Réélu). 


German Region 


The officers for 1958 are: President: O. Heinisch; Secretary-Treas- | 
urer: W. Ludwig. 


Italian Region 


The officers for 1958 are: President: G. Montalenti; Secretary: 
R. Scossiroli; Treasurer: F. Sella. 


Netherlands 


As a result of continual close collaboration with the Medical- 
Biological Section of the Netherlands Statistical Association and the 
Section on Statistical Techniques of the Netherlands Association for 
Agricultural Sciences, members were invited to meetings of these two 
societies. The following papers were read: Precision in clinical medicine 
(Gooszan) ; Precision in clinical chemistry (Strengers); Control systems 
in the clinical laboratory (de Waal; Gooszan); The use of control 
systems (van Strick); Standard control sera (Spaander); Methods of 
increasing precision (van Kampen; Bloemera; Enneking); Probems 
and methods in econometrics (Mol; Liberg) ; Factor analysis and regres- 
sion (Hamming); Tests based on runs (van Strick); Partially balanced 
designs (Justesen); Biometrical genetics (van der Veen); A plant 
breeding problem (Lieben); Errors in soil and plant analyses (Vermeu- 
len); Electronic computers and Monte Carlo methods (Nieden) ; Punch 
cards and computers in medical and biological research (Koene; Muller) ; 
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2 X 2 tables (van Eeden); Statistical regularities in series of observations 
(van Strick); Confidence intervals in industry (Willenze). The group 
heard papers by D. J. Finney on quantal responses and on problems 
of plant selection, and also participated in a discussion on statistical 
research in the human field. 


CHANGES IN MEMBERSHIP 
(FEBRUARY 1958-APRIL 1958) 


Changes of Address 


Mr. James P. Aaron, Jr., 1092 Braddock Road, Cumberland, Marylana, 
U.S. A. 

Mr. Allan Birnbaum, Department of Mathematical Statistics, Columbia 
University, New York 27, N. Y., U.S. A. 

Dr. Charles Bocquet, Laboratoire de Zoologie, Faculte des Sciences, 
Caen, France 

Dr. John W. A. Brant, 3034 Cazador Street, Los Angeles 65, California, 
U.S. A. 

Dr. Lyle D. Calvin, 240 Snell Hall, Oregon State College, Corvallis, 
Oregon, U.S. A. 

Miss Shiela Cohen, ¢/o Rockman, 15 Sokolof, Naharija, Israel 

Dr. Paul M. Densen, Division of Research and Statistics, Health Ins. 
Plan of Greater N. Y., 625 Madison Avenue, New York 3, N. Y., 
U.S. A. 

Mr. Giulio DeVita, Centro Autonomo Militare per l’Energia Nucleare 
(CAMEN), Dipartimento di Radiobiologia, Livorao, Italy 

Mr. Lincoln J. Gerende, 25 Ridgeway, Ann Arbor, Michigan, U. 8. A. 

Mr. James Grizzle, 311 Severin Street, Chapel Hill, North Carolina, 
U.S. A. 

Prof. Dr. Cornelia Harte, Gyrhofstr. 17, Koln-Lindenthal, Germany 

Dr. Robert G. Hoffman, Box 730, J. Hillis Miller Health Center, Uni- 
versity of Florida, Gainesville, Florida, U. 8. A. 

Mr. Alois Kalin, Zederstr. 4, Zurich, Switzerland 

Mr. Jimmy H. K. Kan, Department of Animal Industry, University 
of Arkansas, Fayettesville, Arkansas, U. 8. A. 

Prof. Maxime Lamotte, Laboratoire de Zoologie, Ecole Normale 
Superieur, 24 rue Lhomond, Paris, France 

M. Sully Lederman, Chef de Service a 1’Institut National d’Etudes 
Demographiques, 20-23 Avenue Franklin Roosevelt, Paris, France 

Mr. F. S. McFeely, Department of Mathematics, Montana State 
College, Bozeman, Montana, U. 8. A. 
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Dr. G. G. Meynell, Department of Bacteriology, St. Thomas’ Hospital 
Medical School, London §8.E. 1, England 

Dr. Franklin E. Satterthwaite, 8 Fuller Road, Wellesley 81, Massachu- 
setts, U.S. A. 

Dr. E. F. Schultz, Jr., USDA Southern Research Laboratory, 1100 
Robert E. Lee Blvd., New Orleans, Louisiana, U.S. A. 

Dr. A. R. Sen, 103-A Russa Road, Calcutta-26, India 

Mrs. Ida L. Sherman, 2025 Peachtree Road, Apt. 104, Atlanta, Georgia, 
U.S. A. 

Mr. Raul Vargas, 5514 S. Blackstone Avenue, Chicago 37, Illinois, 
U.S. A. 

Mr. William Weiss, Box 462, R.F.D. 2, Vienna, Virginia, U. S. A. 

Dipl. Ing. Friedrich Wilhelm, Landhaus, Abt. 8, Graz, Austria 


New Members 
At Large 


M. Jose Luis Gonzales, c/o Instituto Mexicano de Investigaciones 
Tecnologicas, del Banco de Mexico, Calzada Legaria 694, Mexico, D.F. 


ENAR 


Mr. Lawrence H. Baker, Department of Animal Husbandry, University 
of Minnesota, St. Paul, Minnesota, U.S. A. 

Mr. Roderick 8. Cowles, Director of Quality Control, E. R. Squibb 
and Sons, Georges Road, New Brunswick, New Jersey, U. 8S. A. 

Mr. William LeRoy Hafley, Institute of Statistics, North Carolina 
State College, Raleigh, North Carolina, U.S. A. 

Mr. Richard A. Lamm, 75 Stewart Manor, Frederick, Maryland, U.S. A. 

Dr. Forrest E. Linder, Director, U.S. National Health Survey, Public 
Health Service, Room 4605 North HEW Building, Washington 25, 
D. C., U.S. A. 

Mr. Laurance G. Locke, 2130 N. Street, NW Apt. 410, Washington, 
D. C., U.S. A. 

Mr. G. B. Oakland, Statistical Laboratory, Department of Agriculture, 
Science Service, Ottawa, Ontario, Canada 

Miss Margaret E. Quinn, 54 Haskell Street, Cambridge 40, Massa- 
chusetts, U.S. A. 


French 


Dr. J. P. Boss, 66 rue des Thermes, Enghien-les-Bains, France 
Mr. P. Corcelle, Institut de Recherches du Coton, 29 rue d’Artois, 
Paris, France 
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Miss A. Gervaise, 2 rue Aunier, Nogent-S/Marne, (Seine), France 

Dr. J. Lejeune, 13 quai Anatole France, Paris 7°, France 

Mr. J. Lellouch, 11 place Adolphe Chérioux, Paris 15°, France 

Mr. H. Malicornet, 11 avenue de Friedland, Paris 8°, France 

Mr. L. Soubies, O.M.I.A., route d’Espagne, Toulouse, France 

Mr. G. Valdeyron, Institut National Agronomique, 16 rue Claude 
Bernard, Paris 5°, France 


German 


Dr. F. Vogel, Max-Planck-Institut fur vergl, Erbbiologie und Erbpath- 
ologie, Berlin Dahlem, Ehrenbergstr. 26-28, Germany 

Prof. Dr. E. Weber, Unter den Linden 6, II Mathematische Institut, 
Humboldt University, Berlin, Germany 


Indian 


Mr. Ram Das, Director, Planning Research and Action Institute, 
U.P., Kalankar House, Lucknow, India 

Prof.S. K. Ekambaram, Head, Department of Mathematical Economics 
and Statistics, University of Mysore, Mysore, India 

Mr. K. P. N. Nair, c/o Department of Research and Statistics, Reserve 
Bank of India, Bombay-1, India 

Mr. R. K. Pande, Statistician, Planning Research and Action Institute, 
Kalankar House, Lucknow, India 

Mr. Narasimha Rao, Assistant Physiologist, Sugarcane Research 
Station, Anakapalle (Andhra Pradesh), India 


Swiss 


Mr. R. Lang, 44 chemin des Fourches, Chene, Bougeries, (Canton de 
Geneve), Switzerland 


WNAR 


Dr. Ole A. Mathisen, Research Associate Professor, Fisheries Research 
Institute, Fisheries No. 5, University of Washington, Seattle 5, 
Washington, U.S. A. 

Dr. Harley Messinger, 3029 Benvenue Avenue, Berkeley 5, California, 
U.S. A. 

Dr. Charles J. Mode, Department of Mathematics, Montana State 
College, Bozeman, Montana, U.S. A. 

Dr. Gordon L. Nordby, P. O. Box 2161, Stanford, California, U.S. A. 
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THE BIOMETRIC SOCIETY 


SECRETARY’S ImMprest ACCOUNT 


Statement of Income and Expenditure during year ending Dec. 31, 1957 


Income £ 8. d. & 
Balance in hand, i Jan. 1957 9 18 #1 
Treasurer—Imprest 535 O 11 
Ottawa Conference fees (1) 3 6 © 

634 19 O 

Expenditure 
Office equipment and stationery 6 6 4 
Secretarial assistance 75 0 0 
Printing 6 6 0 
Postage 37 8 10 125 2 2 

Balance in Hand—31 Dec. 1957 509 17 10 


I certify the above to be a true record of my transactions on behalf of the Bio- 
metric Society. 


24 February, 1958 
M. J. R. Healy (Signed) 
Secretary 


I have examined the account book and vouchers produced by the Secretary 
and certify that the above statement is in avcordance therewith. 


24 February, 1958 
E. Church, A.A.C.C.A. 
(Signed) 


TREASURER’S Report, 1957 
Balance Sheet 


Assets 

Cash: Bank Balance $3,233.20 

Petty Cash 15.50 $3,248.70 

Liabilities 

Subscriptions, 1958 $ 85.75 

Dues, 1958 55.50 $ 141.25 

Surplus $2,539.96 

Gain for period 567.49 3,107.45 


$3,248.70 


Audited: Herbert M. Beckler (Signed) 
Date: April 7, 1958 
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Income and Expenditure Statement 


Income 


Subscriptions—1956 
1957 


Dues —1956 
1957 


Sustaining Memberships—1957 

Back dues and subscriptions 

Consular fee 

Regional allotments 

BIOMETRICS allotments from sustaining members 
Back issues 

Member subscriptions to Journal of A.S.A. 
Overpayments 

Refund from telephone company 

Refund from Post Office 


Less credits and allotments used 
Total Income 


Expenditures 

Postage 

Office supplies and stationery 
BIOMETRICS 

Special services 

Telephone 

Printing 

Salaries 

Second class reentry fee at Knoxville P.O. 
Post Office box rental 

Consular fee 

Secretary’s Office: Directory and operating expenses 
Members subscriptions to Journal of A.S.A. 


Total Expenditures 
Excess of Income over Expenditures 


$ 231.25 


4,520.75 


77.50 


2,098.00 


$ 475.77 
133.35 
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$4,752.00 


2,175.50 


455.00 
28.00 


342.42 


$7,752.92 


$7,185.43 
567.49 


$7,752.92 


Audited: Herbert M. Beckler (Signed) 
Date: April 7, 1958 
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$ 3.00 
a 74.50 
175.00 
68.75 
130.00 
6.70 
2.82 
15.00 
= — 
$ 130.95 
119.70 
5,012.26 
3 46.82 
6.49 
21.61 
184.10 
25.00 
10.50 
3.00 
1,500.00 
125.00 
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OPERATIONS STATEMENT OF BIOMETRICS VOLUME 13 [1957] 


Income (1 February 1958) 
Biometric Society 


566 $ 4.00 $ 2,264.00 
848 2.75 2,332.00 
7 sustaining 25.00 175.00 $ 4,771.00 
547 A.S.A, 4.00 2,188.00 
890 Direct 7.00 6,230.00 
Sale of back issues 
Biometric Society 262.88 
Editor’s Office 4,797.68 5,060.56 
Sale of Reprints 1,338.07 
Exchange 2.64 
$19,590.27 
Expenditures (1 February 1958) 
Cost of Journals 
Printing 
Issue No. 1 $2,930.84 
Issue No. 2 2,313.30 
Issue No. 3 3,152.05 
Issue No. 4 2,586.02 10,982.21 
Mailing and Express Charges 
Issue No. 1 178.66 
Issue No. 2 160.30 
Issue No. 3 211.49 
Issue No. 4 204.89 755.34 11,737.55 
Cost of reprints 
Issue No. 1 * 312.25 
Issue No. 2 i 279.15 
Issue No. 3 280.46 871.86 
Issue No. 4 ———— 
Mailing Charges Reprints 
Issue No. 1 28.17 
Issue No. 2 13.30 
Issue No. 3 12.83 54.30 926.16 
Issue No. 4 ——_- 


Reprinting Six Back Issues 804.12 
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Operating Expenses 

Biometric Society Credits 
Editor’s Office 

Stamps 
Stationery 
Duplication 
Salary and F.I.C.A. 
Office Supplies 
Post Office Fee 
Bank Charges 
Telephone Calls 
Customs 
Addressographing 
Express Charges 


Expenses re. change of Business Office 


Refund on Subscriptions 
Post Office Deposit 


Income 
Expenditure 
Surplus 
Non-operating Items 
Income 
Interest on Bond 
U. 8. Pharmacological Society 
Sale of Society Members List 
Expenditures 
Safe Custody of Bond 
Insurance—back issues 
Surplus 
Gross Surplus 
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Balance Sheet Biometrics Volume 13 


Assets (1 February 1958) 
Accounts Receivable 
Bank Balances 
New York 
Ottawa 
Blacksburg 
U.S. Treasury Bond 
Liabilities 


Subscriptions to Volumes 14, 15, 16 


Balance from previous volumes 
Surplus from Volume 13 


Note: 


$ 21.62 
$ 623.66 
566.58 
57.04 
495.10 
54.34 
10.00 
2.41 
41.98 
28.72 
61.96 
34.89 
492.80 
12.25 
75.25 2,556.98 $ 2,578.60 
16,046.43 
19,590.27 
16,046.43 
$ 3,543.84 
$ 143.75 
200.00 
10.00 $ 353.75 
10.00 
189.70 199.70 
154.05 
$3,697.89 
$ 2,413.96 
$9,079.84 
2,419.28 
2,269.85 13,768.97 
5,000.00 
$ 4,809.50 
12,675.54 
3,697.89 


$21,182.93 $21,182.93 


Not included in assets is stock of back issues from Volumes 1-13 and reprints. 
Not included in expenses is printing bill for December reprints, estimated at 
$300.00 but not received as of 1 February 1958. 


Audited: J. W. Guthridge 
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NEWS AND ANNOUNCEMENTS 


Members are invited to transmit to their National or Regional Secretary 
(if members at large, to the General Secretary) news of appointments, 
distinctions or retirements, and announcements of professional interest. 


L. C. A. Corsten took his doctor’s degree at the Agricultural Uni- 
versity of Wageningen, Holland, on a thesis: vectors, a tool in statistical 
regression theory. 

Dr. C. L. Rumke was admitted as lecturer in statistics at the Kry 
Medical Faculty of the Vrije Universiteit at Amsterdam. On February 
14 he read an introductory paper on the task of the medical statistician. - 


BRUXELLES 1958 
Institut International de Statistique 


Dans le cadre de la 3le Session extraordinaire de |’Institut Inter- 
national de Statistique organisée 4 Bruxelles du 3 au 8 septembre 1958 
se tiendront deux séances conjointes avec la Biometric Society: 

1. Séance organisée par le Professeur G. M. Cox (Raleigh, N.C., 

U.S. A.) Théme: “Recent developments in experimental design.” 

2. Séance organisée par le Dr. L. Martin 

Théme: ‘L’application des méthodes statistiques dans les 
sciences biologiques.” 


Société Adolphe Quctelet 


Journée Biométrique organisée le 9 septembre 1958 par la Société 
Adolphe Quetelet, Région de la Biometric Society pour la Belgique et le 
Congo Belge: 

Thémes: Utilisation des méthodes biométriques, 1) dans les sciences 

agronomiques, 2) dans les sciences médicales, and 3) dans 
les sciences physiques et chimiques. 


Fédération Internationale Pharmaceutique et Société Adolphe Quetelet 


Réunion conjointe de la Fédération Internationale Pharmaceutique 
et de la Société Adolphe Quetelet, le mercredi 10 septembre 1958: 


305 


» 


306 BIOMETRICS, JUNE 1958 


Colloquium sur les applications des méthodes biométriques et 

statistiques aux sciences pharmaceutiques en général. 

Pour tous renseignements complémentaires, s’adresser au Dr. L. 
Martin, 7, rue Héger-Bordet, Bruxelles, Belgique. 


STATISTICAL SUMMER SEMINAR 


The ninth session of the Statistical Summer Seminar will be held 
July 28-August 8, 1958 at the MIT Endicott House, Dedham, Massa- 
chusetts. This seminar is a private, non-profit organization dedicated 
to exploring new areas of application in statistics. It is sponsored by 
the Office of Naval Research. 


The program will be: 


July 28-August 1: “Tolerance Management” with Dr. Emmanuel 

Davis, General Electric Company, as program chairman. 

August 4-8: “Design of Industrial Experiments” with Dr. Martin 

B. Wilk, Bell Telephone Laboratories, as program chairman. 

Each session will have one or two speakers a day as discussion 
leaders. It is felt that this arrangement would provide an unhurried 
atmosphere to foster an exchange and growth of ideas. Familiarity 
through study and work in the above statistical fields or professional 
competance in the field of application is assumed as discussions will 
be at a fairly sophisticated level. 

The registration fee which covers all living expenses is $150 for each 
session. Grants-in-aid are available for educators and government 
workers. Further information may be obtained from Irving Weiss, 
Bell Telephone Laboratories, North Andover, Massachusetts. 


SAMUEL Bropy MEMorRIAL 


A committee has been set up to establish a memorial honoring 
Dr. Samuel Brody who died August 6, 1956. The memorial would be 
a permanent lectureship to be delivered on the University of Missouri 
campus by a distinguished worker in agricultural science. The fre- 
quency of the lectures and their availability in printed form will depend 
on the funds the committee is able to raise. Donations to the memorial 
may be sent to the Samuel Brody Memorial Committee, Eckles Hall, 
University of Missouri, Columbia, Missouri. 
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SyMPosIuM ON BIOMETRICAL GENETICS AND THE 
FourtsH INTERNATIONAL BIOMETRICAL CONGRESS 


The Symposium on Biometrical Genetics and the Fourth Interna- 
tional Biometrical Congress will meet in Ottawa, Canada, from August 
28 to September 2, 1958. Queries regarding programmes and contributed 
papers are to be sent to Mr. M. J. R. Healy, Secretary, The Biometric 
Society, Rothamsted Experimental Station, Harpenden, Herts., Eng- 
land; queries concerning local arrangements should be sent to Mr. G. B. 
Oakland, Statistical Laboratory, Science Service Building, Central 
Experimental Farm, Ottawa, Ont., Canada. A final mailing of informa- 
tion material will be sent to members of The Biometric Society in the 
near future. 
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JOURNAL OF THE 
AMERICAN STATISTICAL ASSOCIATION ,,... 1958 
1757 K St., N. W. Washington 6, D.C. Vol. 53 No. 282 


ARTICLES 

Empiric Investigation of a Test of Homogeneity for Populations 

Composed of Normal Distributions G. A. BAKER 
A New Bivariate Sign Test ISADORE BLUMEN 
Complete Counterbalancing of Immediate Sequential Effects : 

in a Latin Square Design JAMES V. BRADLEY 
Design and Operation of a Double-Limit Variables 

Sampling Plan ACHESON J. DUNCAN 


Effect of Varying Degrees of Transitory Income on Income Elasticity 
of Expenditures MARILYN DUNSING anp MARGARET G. REID 


Factorial Experimentation in Scheffe’s Analysis of Variance 


for Paired Comparisons OTTO DYKSTRA, Jr. 
Actuarial Estimation of Survivorship 

in Chronic Disease LILA ELVEBACK 
A New Index of Manufacturing Production HARRY ERNST 
The Accuracy and Structure of Industry Expectations in Relation 

to Those of Individual Firms ROBERT FERBER 
Graphic Computation of Tau as a Coefficient 

of Disarray HAROLD D. GRIFFIN 
The Precision of Unbiased Ratio-Type 

Estimators H. O. HARTLEY ann LEO A. GOODMAN 
The First 1,945 British 

Steamships J. R. T. HUGHES ann STANLEY REITER 
Inadmissible Samples and Confidence Limits HOWARD L. JONES 
Nonparametric Estimation from 

Incomplete Observations E. L. KAPLAN anp PAUL MEIER 


Weight-Height Standards Based on 
World War II Experience BERNARD D. KARPINOS 


On Noncoverage of Sample Dwellings LESLIE KISH anp IRENE HESS 


The Use of Pea Work ey for Cost Analysis and Control 
A. C. ROSANDER, H. E. GUTERMAN, anp A. J. McKEON 


The Trentile Deviation Method of 
Weather Forecast Evaluation 


BOOK REVIEWS 
PUBLICATIONS RECEIVED 


SUMMARIES OF PAPERS DELIVERED 
AT 117TH ANNUAL MEETING 


THE AMERICAN STATISTICAL ASSOCIATION 
INVITES AS MEMBERS ALL PERSONS INTERESTED IN: 


1. Development of new theory and method 
2. Improvement of basic statistical data 


3. Application of statistical methods to practical problems 


M. J. SLONIN 


For further information, please contact the American Statistical Association, 
1757 K Street, N. W., Washington 6, D. C. 
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INFORMATION FOR CONTRIBUTORS 


MANUSCRIPTS 


Contributions for Biometrics may be addressed to Dr. Ralph A. Bradley, Depart- 
ment of Statistics, Virginia Polytechnic Institute, Blacksburg, Virginia, U.S.A.; 
authors residing in the following Society Regions can expedite consideration of papers 
by submitting them to the appropriate Associate Editor, namely; BRITISH RE- 
GION: Dr. S. C. Pearce, East Malling Research Station, East Malling, Maidstone, 
Kent, England; AUSTRALASIAN REGION: Dr. E. A. Cornish, University of 
Adelaide, Adelaide, Australia; FRENCH REGION: Dr. Georges Teissier, Faculté 
des Sciences de Paris, 1 rue V. Cousin, Paris, France. QUERIES, NOTES, and 
related correspondence should be directed to Professor G. W. Snedecor, Statistical 
Laboratory, Iowa State College, Ames, Iowa, U.S.A. 

MANUSCRIPTS must be submitted in triplicate, with typescript doublespaced 
throughout. Marginal notes may obviate typographical difficulties presented by 
complicated formulae or tables—authors should not attempt editorial instructions 
or markings for the printer. TABLES should be identified by arabic number and 
by a short descriptive title. ILLUSTRATIONS should also be identified by arabic 
number and by a brief caption. (Captions should not be included in illustrations, 
but should be typewritten collectively on an accompanying sheet.) Originals 
should be approximately 8.5 x 11 in. (21.5 x 28 cm.). The original of each chart, 
diagram or graph should be executed in black on white drawing paper or board, on 
blue tracing linen, or on coordinate paper ruled in blue only; coordinate lines to be 
reproduced should be ruled in black. For printing, illustrations may be reduced to 
¥ or \ original dimensions. Lines should therefore be of sufficient thickness, and 
decimal points, periods, and stippled dots should be solid black circles large enough 
to reproduce well. Lettering and numerals should be at least 1 mm. high when 
reproduced in a cut 3 in. (7.5 cm.) wide. Photographs should be prints on glossy 
paper with strong contrasts, and if grouped in a plate should be mounted contig- 
uously. All tables and illustrations should be mentioned explicitly in the text. 
REFERENCES (BIBLIOGRAPHIC) should be collectively listed alphabetically 
by author; textual citation by author and year is preferred. 


ABSTRACTS 


Abstracts of papers presented at meetings of the Biometric Society or of its 
regions are printed in Biometrics following such meetings. They should be submitted 
to the person designated to receive them for a particular meeting in exactly the form 
published in Biometrics (except for an Abstract Number), doublespaced on bond 
paper and in duplicate. Use of formulae requiring display printing is to be avoided. 


Notices, ANNOUNCEMENTS AND Biometric Society Reports 


International and regional reports and notices should be submitted by the 
appropriate officers of the Society and its Regions in duplicate doublespaced on 
separate sheets exactly as they are to be printed in Biometrics. Other material to 
be printed in News and Announcements should also be submitted doublespaced 
and in duplicate. 


Sustarninc MEMBERS OF THE BIOMETRIC Socusrr 


Abbot Laboratories 

American Cancer Society, Inc. — 

Merck, Sharp and Dohme Research Laboratories 
Schering Corporation 

Smith, Kline and French Laboratories 

E. R. Squibb and Sons 

Wyeth Institute of Applied Biochemistry 
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BACK ISSUES 


Back issues of Biometrics are available at the following postage-paid - J 
prices in U.S.A. currency: _— . 


Price per Price per 
Year Volume Number Single Number Volume(unbound) 
1945 1 1 to6 $1.00 $6.00 
1946 2 1 to6 1.00 6.00 
1947 3 lto4 1.50 5.00 
1948 4 lto4 1.50 5.00 
1949 5 1 to4 1.50 5.00 
1950 6 1 to4 1.50 5.00 
1951 7 lto4 2.00 8.00 
1952 8 1 to4 2.00 8.00 
1953 9 lto4 2.00 8.00 
1954 10 lto4 2.00 8.00 
1955 11 2.00 8.00 
1956 12 lto4 2.00 8.00 
1957 13 1 to4 2.00 8.00 


Inquiries, non-member subscriptions, and orders for back issues should 
be addressed to: 

BIoMETRICS 

DEPARTMENT OF STATISTICS 

INSTITUTE 
BuiacksBur@, Virani, U.S.A. 


Reprints of individual articles are not available except +o authors at the 
time of printing. Three special issues are among the numbers listed 
above. They are: 


1947 Volume 3 Number 1 The Analysis of Variance 
1951 Volume 7 Number 1 Components of Variance 
1957 Volume 13 Number 3 The Analysis of Covariance 
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